© L. Orr et al., Published by EDP Sciences 2024

1 Introduction

The European Space Agency (ESA) Swarm mission is a constellation of three polar-orbiting satellites dedicated to measuring the Earth’s magnetic field with unprecedented precision (Friis-Christensen et al., 2008; Hammer et al., 2021). The primary instrumentation consists of a Vector Fluxgate Magnetometer (VFM) to measure the magnetic field strength and direction and an Absolute Scalar Magnetometer (ASM) to measure the absolute magnetic field intensity (Hulot et al., 2015) and to allow cross-calibration of the vector instrument (Toffner-Clausen et al., 2016). The Swarm mission provides a set of daily Level1b products for each satellite containing a time series of corrected, calibrated and quality-assessed measurements, including the scalar and vector magnetic field data at a 1 Hz rate.

The Swarm mission is part of the ESA Earth Explorer series which produces high-quality science data. As such, the established processing chain emphasizes the quality and completeness of the data collection over the speed of delivery. In the case of Swarm, a time lag of four days exists between collection and delivery to the user community. This ensures that all available measurements have been downloaded from the satellites’ onboard memory and that the additional supporting information from other sources required to process the magnetic and electric field data has been captured. Once the Level1b data are released further derived Level2 products can be computed by other institutes (Smith et al., 2022). An example of a daily Level2 product is the British Geological Survey (BGS) ‘fast-track’ magnetospheric model (MMA_SHA_2F) (Hamilton, 2013). The fast-track product is a per-orbit estimate of the degree and order 1 external magnetic field. It is produced daily (hence the ‘fast-track’ moniker), with a four-day lag, providing Gauss coefficients at a 90-min cadence by interpolating the actual per-orbit estimates to a fixed set of timestamps (00:00 UT, 01:30 UT etc.). It is intended to behave similarly to the Disturbance storm time (Dst) Index and characterize the rapidly varying magnetospheric field created by the ring current. The product is computed with the most up-to-date magnetic data and made available within a day of the most recent Level1b data release.

In response to the growing global recognition of the hazard and economic risk posed by space weather (e.g. Oughton et al., 2017; UK Cabinet Office, 2017), ESA has been experimenting with reducing the time taken to release Swarm data to the community. The main obstacles to improving the original four-day lag relate to the availability of ground receiving station bandwidth and overpass opportunities to download satellite data, and the length of the full processing chain to convert data from engineering units to precisely located positions and calibrated SI units. As these bottlnecks have been overcome, the more rapidly released data is known as the ‘FAST’ product and is additional to the standard operational (‘OPER’) Level1b data provided. In theory, it may be possible to provide FAST data on a per-orbit basis with a relatively short time lag of less than an hour. In practice, there is usually a longer wait for the FAST data but consequently, data for several orbits are delivered together.

With the potential of a close-to real-time dataset from across the globe, this study aims to investigate the use of Swarm FAST data for monitoring space weather hazards by designing a Hazard Variation Index (HVI) using 1 Hz magnetic (MAG) data from the three satellites. Such a hazard index would allow maps of localized intense variation of the field to be quickly derived during geomagnetic storms and regions where space weather effects on the ground can be identified. As relatively few regions of the globe have dedicated ground magnetic observatories, even fewer with real-time data capability, our method complements the ground network in areas where in-situ measurements do not exist (e.g. oceans) and measurements are sparse. In such regions, the potential impacts such communications or GNSS loss on long-distance aircraft, marine and space-based infrastructure is of concern.

We propose to establish the range of variation in each magnetic component per gridded cell across the globe. If the standard deviation of the 1 Hz magnetic data in a grid cell exceeds a pre-determined value at the satellite location, this indicates a highly variable ground magnetic field in the region which may produce a potential space weather hazard (i.e. a series of thresholds can be established indicating whether the exceedance implies a low, moderate or high activity). This is similar in theory to determining the deviation away from a standard core field model, or an instantaneous dB/dt variation but is more robust to spikes and offsets, and localises the hazard over a wider region as would be expected for a space weather event. The along-track standard deviation has previously been used to quantify external field signals at satellite altitude and infer ground variability, in order to produce a quantitative, temporo-spatially varying weighting scheme for ground-based observations rather than satellite measurements (Thomson and Lesur, 2007; Brown et al., 2021).

In this study, we describe the calculation of baselines and threshold indices in the Methods section. The Results section provides several examples showing the HVI response to times of low and high geomagnetic activity. The hazard index is shown to highlight localised regions of geomagnetic activity and has a high positive correlation with versions of Kp (Bartels et al., 1939) such as the Hp60 index (Yamazaki et al., 2022); Kp is one of the most commonly used ground-based measures of sub-auroral geomagnetic activity despite its limited time and spatial resolution. We also compare to Hp60 as it is a Kp-like index with a cadence matching closer to that of a half-orbit of the Swarm satellites.

2 Methods

In order to determine whether the standard deviation of the magnetic field measurements at a particular location is anomalous we must first establish the expected variation, at satellite altitude, during relatively normal or quiescent periods. Magnetic field variation has a well-known geomagnetic latitudinal distribution, being smaller at low to mid-latitudes and larger over the auroral and equatorial electrojet zones (Cox et al., 2018). However, the magnitude of the magnetic variation at any given location will have a local time, seasonal and solar cycle dependence (Owens et al., 2021). We use ten years worth of Swarm measurements from 2014–2023 to calculate the background magnetic field variation baselines over almost one solar cycle.

2.1 Data

Swarm magnetic field data files are available from Alpha (A), Bravo (B) and Charlie (C), from the mission started on 25th November 2013, at altitudes of 430–520 km. Swarm Alpha (A) and Charlie (C) fly nearby (side-by-side) at a lower altitude while Swarm Bravo (B) orbits separately at a slightly higher altitude. We use the Swarm OPER 1 Hz Level 1B LR MAG baseline 06 daily files from 01-Jan-2014 to 31-Dec-2023. The orbital planes of A/C and B separate at a rate of about 24°/year (i.e. the difference in longitude of the ascending node) as the satellites precess through all local times at different rates. We will focus on results from Swarm Bravo throughout the manuscript though similar results apply to Swarm Alpha.

The magnetic field data are processed on a per-day basis for convenience. The data were collected via the VirES Python client (Smith et al., 2022) which contains mostly version 0602, with a small selection of 0603 data included. Data with flags above a certain level are excluded from the analysis to avoid including unreliable measurements. Vector data with Flag values of F > 30 (meaning gaps/discrepancies with ASM measurements) and Platform > 67 (meaning gaps/discrepancies Bus and AOCS telemetry), or B = 255 (meaning there are not enough VFM samples to generate B_NEC) or q = 255 (there are not enough star tracker (STR) data to generate attitude information) are excluded from the analysis. Magnetic field intensity is also excluded if F ≥ 16 as this indicates a discrepancy between ASM and VFM measurements.

The VirES system computes magnetic residuals using the CHAOS-7 field model (Finlay et al., 2020) at the satellites’ positions, removing the background field from the measurements of the magnetic field VFM components (B_x, B_y, B_z) and ASM magnetic field intensity (F_s). We remove the core (CHAOS-Core to max degree and order 20), crust (CHAOS-Static to max degree and order 185) and magnetosphere (CHAOS-MMA-Primary [external] and CHAOS-MMA-Secondary [induced] to max degree and order 2) to produce an estimate of the residual ionospheric field (plus unmodelled contributions and noise) as measured by the Swarm satellites.

2.2 Magnetic variation values

We next calculate the 20-s along-track standard deviations to represent the measured magnetic field variations, this is a similar concept to using dB/dt to measure changes in time but more robust to spikes and outliers (Thomson and Lesur, 2007). Note, that we made a comparison of the yearly binned magnetic field variations using 20-s standard deviation with the dB/dt. For Swarm, the two variables show minimal differences in variability hence along-track standard deviation is a suitable parameter for this analysis (see Fig. S1in the Supplementary information).

The magnetic variation values i.e. the along-track standard deviation, ${\sigma }_B({\overline{\theta }}_j,{\overline{\varphi }}_j)$, in each magnetic field component (B represents B_x, B_y, B_z or F_s), are calculated for the mean latitude and longitude of each N = 20 s time window along the satellite orbit:$${\sigma }_B({\overline{\theta }}_j,{\overline{\varphi }}_j)=\sqrt{\sum _{k=1}^N\mathrm{}\frac{{\left(B({\theta }_{i_k},{\varphi }_{i_k})-{\overline{B}}_j\right)}^2}{N-1}},$$(1)

where j is an index of the reduced sample points (daily $j=1,\dots ,\frac{86400}{N}$) and i_k = k+N(j−1) represent the original time points where k∈{1,..,N}. ${\overline{\theta }}_j$, ${\overline{\varphi }}_j$ and ${\overline{B}}_j$ are the mean latitude, longitude and magnetic field at index j, i.e. the mean over the N = 20 time points, i_k.

The choice of time window, N, makes little difference to the global magnetic variation pattern with changes mostly associated with the magnitude. Figure S2 of Supplementary information shows a comparison of the yearly binned magnetic variation with standard deviation taken over N = 5 and N = 10 s, which can be compared to the N = 20-s equivalent in Figure S1. We choose a standard deviation of 20 s as a compromise between computational cost and spatial coverage. 20 s equates to around 150 km, approximately 1 degree of latitude travelled by the satellite.

2.3 Binning using a spherical geodesic grid

Each day’s worth of magnetic variation data, ${\sigma }_B(\overline{\theta },\overline{\varphi })$, is then binned using a spherical geodesic grid based on an icosahedron of triangles (Williamson, 1968; Sadourny et al., 1968). The grid starts as an icosahedron with 20 faces and 12 vertices. A high-performance function (Gagarinov, 2017) then recursively partitions each face into a further 4 triangles with the vertices as the centre of the original face and edges. Each iteration, i, results in 20×4ⁱ faces and 10×4ⁱ+2 vertices. We recursively divide i = 5 times to form a quasi-uniform grid of 20480 faces and 10242 vertices on a unit sphere. Once calculated the coordinates of the triangulation vertices are converted to spherical coordinates and then into latitude and longitude.

Each vertex is the centre of a bin and is indexed with Φ. The bin is the pentagon/hexagon made up of the triangles between the vertex and its neighbouring vertices. This can be imagined as a traditional stitched football, made up of black pentagons and white hexagons (e.g. Alam et al., 2010). Figure 1 shows an example of the gridded vertices as blue dots. The light blue/cyan vertex in the middle of the figure has its associated bin highlighted in orange. The orange triangles between the cyan vertex and its neighbouring orange vertices make up a hexagon. This hexagon is the bin of the cyan vertex and is approximately 1/3405 of the area of a sphere. At a satellite altitude of 500 km with a radius of 6871.2 km, the area is ∼1.75 × 10⁵ km². Figure S3 in Supplementary information shows some further examples of bins associated with various grid points.

Figure 1 Example of standard deviation computation and icosahedron binning. The dark blue dots show the quasi-uniform spherical geodesic grid. The bright blue/cyan dot represents our example grid point with its associated bin shown as the orange hexagon linking it to its neighbour vertices (orange dots). The purple dashed lines represent Swarm satellite passes with purple squares marking the central location of the 20-s standard deviations of the magnetic field. The average of the points (green squares) that fall within this bin are counted. The coastal outline is shown for reference.

2.4 Binned magnetic variation values

When the Swarm satellite passes through the bin the magnetic variations are captured. Over the course of one day the satellite passes that go through the bin are stored and from these a daily mean variation per bin is computed. In Figure 1 the three green squares represent the standard deviation of 20 s of Swarm data (magnetic variation) captured within the bin surrounding the cyan vertex. This includes any ${\sigma }_B(\overline{\theta },\overline{\varphi })$ measured that day within or on the boundaries of the bin. The daily magnetic variation is given by:$${\overline{\sigma }}_B(\mathrm{\Phi },{\psi }_d)=\frac{1}{M}\sum _{j=1}^M\mathrm{}{\sigma }_B({\overline{\theta }}_j,{\overline{\varphi }}_j),$$(2)

where ψ_d is the day of the year, $\{{\overline{\theta }}_j,{\overline{\varphi }}_j\}$ are the latitude and longitude points within bin Φ, and M is the number of values in the bin (e.g. in Fig. 1, M = 4 for 4 green squares). Note the bins of each grid point are overlapping and hence a Swarm magnetic variation measurement ${\sigma }_B({\overline{\theta }}_j,{\overline{\varphi }}_j)$ may be included in the mean of multiple (typically 4–5) vertices ${\overline{\sigma }}_B(\mathrm{\Phi },{\psi }_d)$. The binning results in the daily temporal structure being lost, i.e. we retain information about the mean location but not the associated local time of the measurements.

Additionally the annual magnetic variation is calculated as the mean of the daily variations within a year:$${\overline{\sigma }}_B(\mathrm{\Phi },{\mathrm{\Psi }}_y)=\frac{1}{d_a}\sum _{\psi =1}^{d_a}\mathrm{}{\overline{\sigma }}_B(\mathrm{\Phi },{\psi }_d),$$(3)

where Ψ_y is the year and d_a is the number of days in that year.

Figure 2 shows the annual magnetic variations (Eq. 3) of the magnetic field magnitude determined from Swarm B Level1b 1 Hz data, ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_y)$, for the years Ψ_y = 2014 to 2023 in the spherical geodesic grid (although we show the magnetic field scalar magnitude in Figures 2–5, any component (B_x, B_y, B_y) can be used). While the majority of binned annual mean variation values are ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_y)<0.2$, the variation around the poles is much greater with some bins having average values of ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_y)>2$.

Figure 2 The annual magnetic variations (Eq. 3) of the Swarm B magnetic field intensity residuals.

The increase of the annual magnetic variation compared to quiet time can be calculated with a year at solar minima between Solar Cycle 24 and 25 (taken as Ψ_y = 2020). The normalized annual variation is given by:$$\Delta {\overline{\sigma }}_B(\mathrm{\Phi },{\mathrm{\Psi }}_y)=\frac{{\overline{\sigma }}_B(\mathrm{\Phi },{\mathrm{\Psi }}_y)}{{\overline{\sigma }}_B(\mathrm{\Phi },{\mathrm{\Psi }}_y=2020)}.$$(4)

Figure 3 shows $\Delta {\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_y)$ for the years $\mathrm{\Psi }=$2014 to 2023, where $\Delta {\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_y=2020)=1$ as it has been used as the baseline. For $\Delta {\overline{\sigma }}_{F_s}=1$ (yellow) there is no difference between the yearly variation and that of 2020. For $\Delta {\overline{\sigma }}_{F_s}>1$ (yellow–blue) the yearly variation is high when compared to $\Delta {\overline{\sigma }}_{F_s}=2$ meaning the yearly variation is twice that of 2020. The yearly magnetic variation of 2020 is the lowest overall so, in general, the maps of all years show higher variation ($\Delta {\overline{\sigma }}_{F_s}>1$). The years 2015, 2022 and 2023 are the most active with $\Delta {\overline{\sigma }}_{F_s}>3$ in places i.e. the yearly variation is more than three times that of 2020. On average, 2023 has a high variation compared to the previous nine years due to the approach of the peak of the solar cycle. Figure S4, Supplementary information shows the same as Figure 3 but uses 2023 to normalize against, as the year with the greatest magnetic variation.

Figure 3 The normalized annual magnetic variation of ${\overline{\sigma }}_{F_s}$ according to equation (4), with that of the quiet year 2020, using magnetic field intensity measurements from Swarm B as in Figure 2. $\Delta {\overline{\sigma }}_{F_s}=1$ (yellow) when the yearly binned variation is equal to that of 2020. $\Delta {\overline{\sigma }}_{F_s}>1$ (green–blue) when the variation is greater than 2020 and $\Delta {\overline{\sigma }}_{F_s}<1$ (orange) when the variation is less than in 2020.

2.5 Quantiles determined from ten years of daily magnetic variation

We compute the quantiles of daily standard deviation based on ten years of Swarm data. The exceedance quantile values, q_n,B, of the daily variations of the magnetic field components, ${\overline{\sigma }}_B(\mathrm{\Phi },{\psi }_d)$ are calculated for n = 50, 60, 70, 80, 90, 95, 97, 98 and 99, for the years 2014–2023. To do this, the ten years of daily magnetic variation values were concatenated, sorted and the nth percentile per bin was found. For Swarm B the bins contain a minimum of 654 data points, whereas A has a minimum of 788 per bin. A and B respectively have a maximum of 3634 and 3642 data points per bin (there are 3652 days in ten years), with an average of 1403 for A and 1148 for B. As a very rough approximation, if the bins span 4° longitude and there are 15 full orbits in a day we would expect approximately a third of the bins to pass through daily ($\frac{4\times 30}{360}=\frac{1}{3}$), with a lot of variation depending on latitude.

Figure 4 shows q_n,Fs(Φ), the quantile values in each bin, with a narrow colour range to highlight spatial variation in the lower percentiles. Higher values are seen at the poles and equator. With the quantile maps, the exceedance of the nth quantile per component can be directly looked up and used as thresholds to indicate whether the instantaneous magnetic variation value (e.g. of ${\sigma }_B(\overline{\theta },\overline{\varphi })$ from Eq. 1) implies low, moderate or high geomagnetic activity at that location.

Figure 4 The exceedance quantile values, $q_{n,F_s}$, of the daily magnetic variations, ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\psi }_d)$ for the years 2014–2023 from the Swarm B magnetic field intensity residuals. The ten years of daily magnetic variation values were concatenated, sorted and the nth percentile per bin is found.

As an approximation, to work out how often low, moderate or high geomagnetic activity occurs, a calculation of the statistics of the Kp index from 1st January 1932–2nd April 2023 using data from Matzka et al. (2021) is shown in Table 1. Large storms (Kp≥6) occur approximately 1.68% of the time.

A selection of rank order plots is shown in Figure S5 of the Supplementary information to show the distributions of the daily variation, ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\psi }_d)$ and the quantile values q_n,Fs at various latitude bins.

For reference, Table 2 shows the equivalent Kp value for the quantiles chosen for Figure 4. Even at the 99% percentile, Kp is just below 7. We note that Chakraborty and Morley (2020), who studied the statistics of Kp, also chose their threshold for storm time as Kp≥5⁻, which is equivalent to our 95% percentile.

Using the quantile percentages in Table 2 as a reference for high or low activity in terms of our defined quantiles, we assume ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },\psi )\ge q_{95,F_s}\left(\mathrm{\Phi }\right)$ is indicative of Kp≥5⁻ and ‘storm time’.

2.6 Defining an exceedance parameter

With the established quantile thresholds we can then calculate the fraction of bins where the daily magnetic variation, ${\overline{\sigma }}_B(\mathrm{\Phi },{\psi }_d)$ exceeds each nth quantile. We name this exceedance parameter Ω_n,B, which is defined as:$${\mathrm{\Omega }}_{n,B}\left(\psi \right)=\frac{1}{m}\sum _{i=1}^m\mathrm{}{\omega }_{n,B}(i,\psi ),$$(5)

where m is the number of bins and,$$\begin{array}{lll}{\omega }_{n,B}(\mathrm{\Phi },{\psi }_d)& =& 1,\hspace{1em}\mathrm{if}\hspace{1em}{\overline{\sigma }}_B(\mathrm{\Phi },{\psi }_d)\ge q_{n,B}\left(\mathrm{\Phi }\right),\\ {\omega }_{n,B}(\mathrm{\Phi },{\psi }_d)& =& 0,\hspace{1em}\mathrm{otherwise}.\end{array}$$(6)

Further, we calculate a measure of the mean exceedance parameter which we will define as the quantile index, QI_B(ψ). It is intended to summarise geomagnetic activity over a defined period with a single parameter and is defined as:$$Q{I}_B\left(\psi \right)=\sum _n\mathrm{}{\mathrm{\Omega }}_{n,B}\left(\psi \right),$$(7)

where n = [50, 60, 70, 80, 90, 95, 97, 98, 99]. The quantile index, QI_B(ψ) is a dimensionless parameter between 0 and 9 per day (ψ), indicating how many of the n thresholds have been exceeded. Equations (5)–(7) can also be calculated per hour or per orbit, as will be shown below. Starting with the lower panel of Figure 5 shows the daily exceedance parameter, ${\mathrm{\Omega }}_{n,F_s}$, associated with the magnetic field intensity for the year 2017. By design, ${\mathrm{\Omega }}_{0,F_s}=1$ i.e. all binned data are equal to or above the 0th quantile. Over-plotted on top of ${\mathrm{\Omega }}_{0,F_s}$, is ${\mathrm{\Omega }}_{50,F_s}$ (lightest purple colour) representing the fraction of bins with a daily value above the 50th quantile, $q_{50,F_s}$ e.g. 80% of bins or 0.8 as a fraction. The shortest daily bars (darkest green) represent the ${\mathrm{\Omega }}_{99,F_s}$ fraction of bins with a daily value above the 99th quantile, $q_{99,F_s}$ e.g. 5% or 0.05.

Figure 5 Fs indices for the year 2017 using Swarm B Level1b OPER data. The top panel overplots the daily quantile index, QIB(ψ) (Eq. 7) with C9 (a daily estimate of Kp). The bottom panel shows the daily exceedance parameter, ${\mathrm{\Omega }}_{n,F_s}$, associated with the daily magnetic field intensity for the year 2017 (Eq. 5). Overplotted coloured bars indicate the proportion of bins in exceedance of the nth quantile. See text for further details.

Moving to the upper panel of Figure 5 the orange line is the daily C9 index shown for reference. C9 is a qualitative estimate of the overall level of magnetic activity for the day determined from the sum of the daily ap amplitudes and converted to be between 0 and 9, similar to Kp (Dieminger et al., 1996). The panel demonstrates that the daily quantile index $Q{I}_{F_S}$ tracks and correlates strongly to C9.

The height of the $Q{I}_{F_S}$ bars in the upper panel tracks the variation of the light-coloured bars in the lower panel, generally corresponding to the C9 values. For example, the largest storm of the year, the 7–8th September 2017 shows a high proportion of bins exceeding their 99th threshold. The effects of heightened solar activity and recurrent coronal holes can be seen from the approximately 27-day periodicity.

3 Application to geomagnetic storms

We next examine the methodology with Swarm data as measured during some of the largest geomagnetic storms of the mission lifetime to date. The quantile threshold exceedance can be applied to either the operational (OPER) dataset or the fast track (FAST) Swarm Level1b product files. To do this we return to the magnetic variation parameter, ${\sigma }_B(\overline{\theta },\overline{\varphi })$ (Eq. 1). Here we show only the scalar component of the magnetic field F_s but B_x, B_y, B_z and F_v are shown in Supplementary information, Figures S6–S9.

3.1 7–8th September 2017

Although the geomagnetic storm of the 7/8th September 2017 was moderate with a Dst of ∼−150 nT, it was still one of the largest events in Solar Cycle 24 producing observed space weather impacts at mid to high latitude (e.g. Dimmock et al., 2019). Figure 6(a) plots ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })$ over the time period from 07-Sep-2017 18:00 UT to 09-Sep-2017 06:00 UT for all of the Swarm B orbits during the storm. The coloured dots represent the magnitude of the magnetic variation and the location of the markers are the mean latitude and longitude of the 20-s standard deviation time windows.

Figure 6 Hazard variation plots for the first peak of the 7–8th September 2017 storm using Swarm B Level1b OPER data. (a) plots the magnetic field intensity variation ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })$ over the time period. (b) shows $\Delta {\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_{2020})$ i.e. panel (a) normalized by the mean magnetic variation of 2020. (c) plots the exceedance of the quantiles $q_{n,F_s}(\overline{\theta },\overline{\varphi })$ where the colour scheme reflects if ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })\ge q_{n,F_s}(\overline{\theta },\overline{\varphi })$. (d) and (e) are as in Figure 5, except each timepoint represents when Swarm B passes between ±75° latitude and (d) over plots Hp60 to better match the cadence of these orbit segments. (f) and (g) are of the form of (c) but for the components Bx, By and Bz, respectively. They only show the exceedance of the highest quantiles (n ≥ 95).

Figure 6(b) shows $\Delta {\overline{\sigma }}_B(\mathrm{\Phi },{\mathrm{\Psi }}_y)$ i.e. the magnetic variation ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })$ values during the storm normalized by ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_y=2020)$, the quietest yearly magnetic variation values (Eq. 4). However, ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_y=2020)$ is the binned yearly magnetic variation and to directly compare with ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })$ it is necessary to interpolate the bin location (Φ) back into the mean latitude and longitude ($\overline{\theta }$ and $\overline{\varphi }$) measured by the Swarm satellite during this time period. The triangulation-based linear interpolation method fits a surface of the form $\overline{\theta }$ and $\overline{\varphi }$ to the binned data. The difference with the quiet year is then defined as:$$\Delta {\sigma }_{F_s}(\overline{\theta },\overline{\varphi })=\frac{{\sigma }_{F_s}(\overline{\theta },\overline{\varphi })}{{\overline{\sigma }}_{F_s}(\overline{\theta },\overline{\varphi },{\mathrm{\Psi }}_y=2020)}.$$(8)

The year 2020 was the lowest activity year and hence, panel 6(b) shows whether ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })$ is higher or lower in comparison. Many of the orbits show the magnetic variation to be more than three times that of 2020. Figure 6(c) plots the exceedance of the quantiles $q_{n,F_s}(\overline{\theta },\overline{\varphi })$ for n = [50, 60, 70, 80, 90, 95, 97, 98, 99] where the colour scheme reflects if ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })\ge q_{n,F_s}(\overline{\theta },\overline{\varphi })$. As in panel 6(b), the binned quantiles ($q_{n,F_s}\left(\mathrm{\Phi }\right)$) have been interpolated to latitude and longitude ($q_{n,F_s}(\overline{\theta },\overline{\varphi })$).

Figure 6(d) plots the HP60, an hourly version of the Kp index (Yamazaki et al., 2022), on the left axis while on the right-hand axis the quantile index per half-orbit, $Q{I}_{F_s}\left(\alpha \right)$. This is calculated from equations (5)–(7) per half-orbit (defined as pole to pole) instead of the day (ϕ=α) and for the mean latitude and longitude instead of the bin ($\mathrm{\Phi }=[\overline{\theta },\overline{\varphi }]$), using the interpolated quantiles $q_{n,F_s}(\overline{\theta },\overline{\varphi })$. Note for 6(d, e) we only include latitude values of $\left|\overline{\theta }\right|<75$ to exclude the response at the poles.

Figure 6(e) is the exceedance parameter per orbit, ${\mathrm{\Omega }}_{n,F_s}\left(\alpha \right)$, where the magnetic variation of F_s components, ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })$ is above the nth quantile, $q_{n,F_s}(\overline{\theta },\overline{\varphi })$ (Eqs. 5 and 6). On the bottom row, panels 6(f–h) are equivalent to the exceedances plotted in 6c but for B_x, B_y, B_z, and only n = [95, 97, 98, 99], i.e. the high quantiles representing approximately ‘storm time’ with Kp ≥ 5⁻ (from Table 2). Note the plots show the B_x, B_y, B_z 20-s standard deviations of the residuals which are not directly proportional to F_s.

The time frame chosen for Figure 6 highlights a period of relative quiet before the storm plus the two peaks and the beginning of the recovery phase. This allows us to see both the quiet orbits and the active ones. According to the description by Walach and Grocott (2019), the Sym-H index shows the main phase of the storm began at 23:07 UT on 7 September 2017. We see an increase in measured activity from ∼22:30 UT with maximum per-orbit variations observed around ∼23:30 UT. During the second peak, we see a very sharp increase in variation around ∼12:00 UT. The maps from 6(a–c) show that the external magnetic field variation varies by orbit, with much higher variation in the west and around 0° longitude. The magnetic field vector components (B_x, B_y and B_z) in 6(f–h) show few exceedance values at low latitudes and again for longitudes 50–180°. High values suggest the auroral region experienced the largest magnetic field variation. In the vector components, all components see a similar number of grid points exceeding the 95th percentile.

3.2 Application with the Swarm FAST data product

In Figure 7, we have used data from the Swarm Level1b FAST product instead of the OPER version. The FAST data are intended to be available as soon as possible for a number of previous orbits which varies depending on how long it has been since the satellite last provided data to its downlink. In general, FAST and OPER magnetic values are identical but the length and temporal coverage of the FAST product varies from file to file. Figure S10, Supplementary information, shows the same figure but using the OPER data instead of FAST. For Figure 7, we show an example from 5th November 2023, a recent storm that had a peak for Kp 7+. The panels have the same layout as Figure 6.

Figure 7 Hazard variation plots for the 5th November 2023 using Swarm B Level1b FAST data. Figure is of the same format as Figure 6.

In panels (a–c), ${\sigma }_{F_s}$ values begin to increase from 10:00 UT with a peak exceedance between 16:00 and 20:00 at longitudes of both ∼150° and ∼−30°, with about 80% for grid points between ±75° latitude exceeding the 95th percentile (our threshold for ‘storm time’) and <60% the 99th. Panel 7(f) shows particularly strong variation in the B_x component in the west.

The quantile index, $Q{I}_{F_s}\left(\alpha \right)$, tracks Hp60 throughout the period, starting at $Q{I}_{F_s}\sim 1$ when Hp60 = 2, increasing the $Q{I}_{F_s}\sim 7$ when Hp60 = 7 before gradually decreasing through the remainder of the time period. $Q{I}_{F_s}\left(\alpha \right)$ then gradually decrease from this peak at ∼17:30 UT for the remainder of the period, roughly following the return to lower HP60. Further examples using FAST Swarm B with low-moderate activity from 1st–5th February 2024 can be found in the Supplementary information, Figures S11–S15.

4 Discussion and conclusions

We have developed and tested a method with Swarm OPER data for use with Swarm FAST data (or indeed any other low Earth orbit satellite mission) by designing a Hazard Variation Index (HVI) using the Level1b 1 Hz LR MAG values. When the value of the standard deviation in any magnetic component exceeds a series of predetermined thresholds, this indicates a highly variable magnetic field in the region which may produce a potential space weather hazard on the ground, not just at orbit altitude.

While this is a similar idea to dB/dt variation or variation away from a core field model, it is more robust to spikes and localizes the hazard over a wider region as would be expected for a space weather event. It also contextualizes the variation with respect to the expected level of magnetic activity. As more satellite data become available, spanning over one solar cycle, the background percentiles will become more robust and comparison between solar cycles can be analysed too.

We point out that the loss of local time separation does affect the results and is particularly clear on the day vs night side of the planet as in Figure 6(d) and Figure 7(d) where the orange line, $Q{I}_{F_s}$, see-saws up and down from orbit to orbit. It would be possible but a much greater effort to sort the input datasets into local time, to determine the HVI. This would likely remove the see-saw effect of the day versus nightside of the planet but is unlikely to change the overall behaviour of the HVI significantly.

To apply the technique to other satellite missions, it is likely each will have to be ‘trained’ separately as the noise floor of each instrument suite may be different. As Swarm is a very high-quality mission, it sets the baseline level for noise. Cube-sat-type missions or platform magnetometers could also be used once their noise characteristics are determined (Olsen et al., 2020).

In terms of application, space weather hazards can have impacts on aircraft HF communications, GNSS accuracy and operation of low Earth orbit satellites. Although true real-time measurements from magnetic satellites (e.g. within 5 min) are not presently available, delivery of Swarm FAST data within an hour and computation of the HVI would be useful to provide post-facto analysis for remote regions. Any technology experiencing unusual effects in polar or oceanic areas could quickly determine if local enhanced space weather activity was a factor. While current coverage only includes the three Swarm satellites, more satellites could be included (e.g. MSS-1a, NanoMagSat). This study is a proof of concept that could be expanded with further magnetic satellite missions, yet still provides better spatial coverage than anything possible currently with ground-based data in terms of low-spatial-resolution indices that are commonly used to infer geomagnetic activity.

We have shown the HVI tracks geomagnetic activity (as measured by the Kp or other global indices) both in quiet and active periods using both OPER and FAST Swarm magnetic data products. The method is quite general and can be used with any good-quality magnetic satellite data feed.

Acknowledgments

We thank ESA for providing Swarm data to the community. Access to Swarm data is via: https://swarm-diss.eo.esa.int/. We use the Python package, viresclient (Smith et al., 2024b), to access OPER and FAST Level1b 1 Hz LR MAG from ESA’s VirES for Swarm service (Smith et al., 2024a). The results presented in this paper rely on data collected at magnetic observatories. We thank the national institutes that support them and INTERMAGNET for promoting high standards of magnetic observatory practice (www.intermagnet.org). GFZ provided Matlab API for download of Kp and derived indices, available via: https://kp.gfz-potsdam.de/en/data (Matzka et al., 2021). This paper is approved by the Director of the British Geological Survey (UKRI).

Supplementary information

Figure S1: Swarm B 2014 yearly mean comparing mean binned standard deviation with mean binned dB/dt over 20 s intervals.

Figure S2: Swarm B 2014 yearly mean comparing mean binned standard deviation over 5 and 10 s intervals.

Figure S3: Two further examples the icosahedron binning in the same style as Figure 1, main text. The left is a vertex with only 5 neighbours forming a pentagon bin (over the Red Sea for reference). The right shows two separate satellite passes being included in the bin (over South England for reference).

Figure S4: Swarm B yearly Fs normalized by the active year of 2023 in the same format as Figure 3, main text.

Figure S5: Swarm B distribution of ${\sigma }_{F_s}$ across the years 2014–2023.

Figure S6: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the Bx component.

Figure S7: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the By component.

Figure S8: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the Bz component.

Figure S9: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the Bv component.

Figure S10: Hazard variation plots for the 5th November 2023 using OPER Swarm B measurements, in the same format as Figure 7, main text.

Figure S11: Hazard variation plots for the 1st February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.

Figure S12: Hazard variation plots for the 2nd February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.

Figure S13: Hazard variation plots for the 3rd February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.

Figure S14: Hazard variation plots for the 4th February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.

Figure S15: Hazard variation plots for the 5th February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.

References

All Tables

All Figures

Figure 1 Example of standard deviation computation and icosahedron binning. The dark blue dots show the quasi-uniform spherical geodesic grid. The bright blue/cyan dot represents our example grid point with its associated bin shown as the orange hexagon linking it to its neighbour vertices (orange dots). The purple dashed lines represent Swarm satellite passes with purple squares marking the central location of the 20-s standard deviations of the magnetic field. The average of the points (green squares) that fall within this bin are counted. The coastal outline is shown for reference.

In the text

	Figure 2 The annual magnetic variations (Eq. 3) of the Swarm B magnetic field intensity residuals.
In the text

Figure 3 The normalized annual magnetic variation of ${\overline{\sigma }}_{F_s}$ according to equation (4), with that of the quiet year 2020, using magnetic field intensity measurements from Swarm B as in Figure 2. $\Delta {\overline{\sigma }}_{F_s}=1$ (yellow) when the yearly binned variation is equal to that of 2020. $\Delta {\overline{\sigma }}_{F_s}>1$ (green–blue) when the variation is greater than 2020 and $\Delta {\overline{\sigma }}_{F_s}<1$ (orange) when the variation is less than in 2020.

In the text

	Figure 4 The exceedance quantile values, $q_{n,F_s}$, of the daily magnetic variations, ${\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\psi }_d)$ for the years 2014–2023 from the Swarm B magnetic field intensity residuals. The ten years of daily magnetic variation values were concatenated, sorted and the nth percentile per bin is found.
In the text

Figure 5 Fs indices for the year 2017 using Swarm B Level1b OPER data. The top panel overplots the daily quantile index, QIB(ψ) (Eq. 7) with C9 (a daily estimate of Kp). The bottom panel shows the daily exceedance parameter, ${\mathrm{\Omega }}_{n,F_s}$, associated with the daily magnetic field intensity for the year 2017 (Eq. 5). Overplotted coloured bars indicate the proportion of bins in exceedance of the nth quantile. See text for further details.

In the text

Figure 6 Hazard variation plots for the first peak of the 7–8th September 2017 storm using Swarm B Level1b OPER data. (a) plots the magnetic field intensity variation ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })$ over the time period. (b) shows $\Delta {\overline{\sigma }}_{F_s}(\mathrm{\Phi },{\mathrm{\Psi }}_{2020})$ i.e. panel (a) normalized by the mean magnetic variation of 2020. (c) plots the exceedance of the quantiles $q_{n,F_s}(\overline{\theta },\overline{\varphi })$ where the colour scheme reflects if ${\sigma }_{F_s}(\overline{\theta },\overline{\varphi })\ge q_{n,F_s}(\overline{\theta },\overline{\varphi })$. (d) and (e) are as in Figure 5, except each timepoint represents when Swarm B passes between ±75° latitude and (d) over plots Hp60 to better match the cadence of these orbit segments. (f) and (g) are of the form of (c) but for the components Bx, By and Bz, respectively. They only show the exceedance of the highest quantiles (n ≥ 95).

In the text

	Figure 7 Hazard variation plots for the 5th November 2023 using Swarm B Level1b FAST data. Figure is of the same format as Figure 6.
In the text

	Figure S1: Swarm B 2014 yearly mean comparing mean binned standard deviation with mean binned dB/dt over 20 s intervals.
In the text

	Figure S2: Swarm B 2014 yearly mean comparing mean binned standard deviation over 5 and 10 s intervals.
In the text

	Figure S3: Two further examples the icosahedron binning in the same style as Figure 1, main text. The left is a vertex with only 5 neighbours forming a pentagon bin (over the Red Sea for reference). The right shows two separate satellite passes being included in the bin (over South England for reference).
In the text

	Figure S4: Swarm B yearly Fs normalized by the active year of 2023 in the same format as Figure 3, main text.
In the text

	Figure S5: Swarm B distribution of ${\sigma }_{F_s}$ across the years 2014–2023.
In the text

	Figure S6: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the Bx component.
In the text

	Figure S7: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the By component.
In the text

	Figure S8: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the Bz component.
In the text

	Figure S9: Hazard variation plots for the 7–8th September 2017 using OPER Swarm B measurements, in the same format as Figure 7, main text, but focusing on the Bv component.
In the text

	Figure S10: Hazard variation plots for the 5th November 2023 using OPER Swarm B measurements, in the same format as Figure 7, main text.
In the text

	Figure S11: Hazard variation plots for the 1st February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.
In the text

	Figure S12: Hazard variation plots for the 2nd February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.
In the text

	Figure S13: Hazard variation plots for the 3rd February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.
In the text

	Figure S14: Hazard variation plots for the 4th February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.
In the text

	Figure S15: Hazard variation plots for the 5th February 2024 using fast Swarm B measurements, in the same format as Figure 6, main text.
In the text