© V. Lanabere et al., Published by EDP Sciences 2024

1 Introduction

Extreme space weather events that represent a real threat to ground- and space-based infrastructure are categorized as low-frequency but high-consequence events, resulting in significant social and economic impacts (Eastwood et al., 2017). These extreme events can have consequences due to a multitude of factors that vary depending on the nature of the impact. For example, some key factors including the location of the system (space/ground-based), physical properties of the hardware, time of the events, and the nature in which the external driver couples to near-Earth space are all important. For instance, during a strong solar event on 3–5 April 2010 the Intelsat satellite G-15 stopped responding to ground commands (Allen, 2010; Loto’aniu et al., 2015), although it was able to reboot itself and restore control later that year. Similarly, on 3 February 2022, SpaceX lost 38 out of 49 Starlink satellites due to increased neutral atmospheric density associated with a geomagnetic storm (Fang et al., 2022; Kataoka et al., 2022). However, there are clear extreme events that produce negative effects in almost all the infrastructures mentioned. For instance, in the 21st century the Halloween storm, peaked between 29 and 30 October 2003, with a minimum Dst of about −400 nT. This storm caused multiple satellites to be affected, including a Japanese satellite that was lost completely; severe high-frequency radio blackout that affected commercial airlines; excessive radiation exposure for high-altitude aircraft and astronauts in the International Space Station (Webb and Allen, 2004; Doherty et al., 2004; Xue et al., 2023). At ground level, geomagnetically induced currents (GICs) in southern Sweden caused a regional blackout that affected the city of Malmö (Pulkkinen et al., 2005). Additionally, airlines rerouted polar flights due to poor high frequency communications (Webb and Allen, 2004). A storm of this magnitude is estimated to occur approximately once every 50 years (Bergin et al., 2023). Recently, from 10 to 12 May 2024, a strong geomagnetic storm with a Dst value of about −412 nT was observed. However, there are still no official reports on how the storm affected different infrastructures. Nevertheless, more significant events were observed during the 20th century, such as those in March 1989 (Dst = −589 nT) and May 1921 (Dst −900 nT) (Kappenman, 2006; Hapgood, 2019). The largest known geomagnetic storm occurred during the 19th century known as the Carrington event, with an estimated Dst of approximately 900 nT (Cliver and Dietrich, 2013).

The potential risk of disruption to ground-based infrastructure during extreme geomagnetic storms has prompted several studies to identify regions where the geoelectric fields are strong enough to generate hazardous GIC. In a recent study by Love et al. (2022), the authors performed a risk analysis of the United States power systems during an extreme event on 13–14 March 1989, which caused a blackout in Québec, Canada, and damaged a high-voltage transformer in the northeastern United States. The authors constructed geoelectric field maps and compared the field’s magnitude with reported power system anomalies, emphasizing the important role of surface impedance in affecting the amplitudes of geoelectric fields. Similarly, Malone-Leigh et al. (2024) generated geoelectric field hazard maps across Ireland using magnetotelluric transfer functions to model the geoelectric field over a 28-year period using ground-based magnetometer data. By estimating GIC in Portugal, Alves Ribeiro et al. (2021, 2023) identified the critical substations in the Portuguese power network that are prone to the strongest GIC effects. Such studies offer valuable insights to support the planning and design of more resilient transmission systems.

Recordings of extreme events in different regions of the geospace environment are generally not available. However, extreme value (EV) theory provides a framework for estimating unobserved levels using historical data (Coles, 2001). The extreme value distributions commonly used are the generalized extreme value (GEV) and Generalized Pareto (GP) distributions, which are specified by three parameters: location (u), scale (σ), and shape (ξ). The ξ parameter characterizes the heaviness of the distribution’s tail; larger values indicate heavier tails, implying that extreme events are more likely to occur. In the field of space physics, extreme value analysis of energetic particles in the van Allen radiation belt (Koons, 2001; Meredith et al., 2015, 2016, 2017) can be useful for spacecraft engineering and evaluating realistic disaster scenarios. Extreme value analysis of geomagnetic indices has also been conducted by several authors. For example, Koons (2001) was the first to apply EV theory to the Ap index, and later, Tsubouchi and Omura (2007) applied it to the Dst index. More recently, studies using new indices and longer databases (Love, 2021; Bergin et al., 2023) have been useful for risk assessment and planning mitigation strategies against space weather hazards. In particular, expected return values of E were used by Patterson et al. (2023) to study the impact of GICs on the electrified railway signaling system in the United Kingdom.

Geomagnetic storm-time GICs have been observed to affect power transmission systems located at various latitudes (Trivedi et al., 2007; Gaunt and Coetzee, 2007; Watari et al., 2009; Pulkkinen et al., 2005). The magnitude of the time derivative of the geomagnetic field |dB/dt| is often used as a proxy to study GICs in power systems. Due to the ground’s electrical conductivity properties and Faraday’s law of induction, an induced geoelectric field is settled, causing currents to flow in long conducting systems such as power grids (e.g. Pulkkinen et al., 2017). Recognizing the importance of 65 and 75° national preparedness to GICs, several EV analyses of |dB/dt| and |E| quantities have been conducted. As an example, Thomson et al. (2011) utilized 1-min ground geomagnetic field data from 28 stations in various locations across Europe, including just two stations in Sweden. The authors employed the peak-over-threshold method and declustered data by identifying extreme values that were separated by at least 12 h. The results showed that stations within geomagnetic latitudes 53–62°N can experience higher extreme magnitudes compared to those predicted between 65 and 75°N, and the shape parameter is positive but generally below one. A similar study in Europe was later conducted by Wintoft et al. (2016). In this study, 1-min geomagnetic data from 14 stations, including two in Sweden, were analyzed using EV analysis applied to dB/dt and the 1-D estimated E. The authors selected extreme events by taking the maximum values from one-year time blocks. Their results indicate that stations south of 60°N magnetic latitude have positive shape parameters, while north of 60°N have shape parameters close to zero. Another EV analysis was conducted in Ireland by Fogg et al. (2023). The study examined the expected 1 in 50- and 100-years magnitude of B and dB/dt and emphasized the importance of more localized space weather predictions at a national level. These earlier studies on EV analysis of dB/dt and E have shown to be very important to understanding the levels and frequency of extreme events.

The Swedish electrical power transmission system is occasionally affected by high-amplitude GICs due to its high geomagnetic latitude location (55°–66.5°N), where the largest geomagnetic variations are recorded (Schillings et al., 2022). This aligns with the results of previous studies, which found larger magnetic perturbations to be expected at Sweden’s latitude. However, these studies have only utilized a limited number of stations in Sweden and/or low-resolution data. These two factors should be taken into account when conducting a study in Sweden. First, it is important to include more stations; this is particularly important since it was shown that the magnitude of the time derivative of the horizontal magnetic field |dH/dt| varies significantly over spatial scales of hundreds of kilometers. Even when |dH/dt| is high over a given region, individual stations within that region can measure significantly different values (Dimmock et al., 2020). Moreover, ground conductivity in Sweden varies by an order of magnitude, and the E thresholds will be dependent on that. Second, calculating the geoelectric field based on 1-min data may result in an underestimation of its peak value to some extent. For instance, Pulkkinen et al. (2006) investigated the degradation of GIC estimates with decreasing sampling rates during the geomagnetic storm period of October 29–31, 2003. The authors observed that the correlation between the modeled geoelectric fields at a 1-s and 60-s sampling rate was approximately 0.9. However, a more recent study by Gannon et al. (2017) examined five geomagnetic storm periods and determined that shorter period magnetic field variations (ranging from 10 s to 1 min) could lead to a 50% loss in peak value compared to 10-s data. Thus, these previous studies show that the degradation of the E estimation depends on the event and the location. In this study, we performed an EV analysis of E in Sweden using 10-s data and a high-resolution conductivity model. Next, we calculated the maximum expected values in 10-, 50-, and 100 years at different locations across Sweden, which are important for national preparedness for GICs.

The purpose of this study is to conduct an EV analysis on high-resolution |E| data in Sweden at 18 locations distributed across six different latitudes to address the following questions: 1) what is the best way to construct independent extreme data points?, 2) is there any difference between different declustering methods?, 3) is the distribution of extreme events the same at different latitudes in Sweden?, 4) what is the maximum expected |E| in 10-, 50-, and 100 years at different latitudes in Sweden?, 5) does the maximum |E| during past events correspond to a 1 into 10-, 50-, or 100-years event? The manuscript is structured as follows. In Section 2, we describe the data, methods, and model used to compute the E in Sweden. In Section 3, we present two different declustering methods and explore different methods for selecting a high threshold to construct independent extreme events. In Section 4, we fit the Generalized Pareto distribution function to two different extreme data sets, compare the obtained parameters, and compute the maximum expected |E| in 10-, 50-, and 100-years. The results are compared with Swedish stations in Section 5. A discussion of the results and future implications are presented in Sections 6 and 7, respectively. Finally, a summary and conclusion are given in Section 8.

2 Data and models

A long-time series of geoelectric field data in Sweden requires a correspondingly long time series of the ground magnetic field data from the region. However, only three stations located in Sweden provide more than 11 years of data. Therefore, the 2-D Spherical Elementary Current System (SECS) method (Vanhamäki and Juusola, 2020) which is used for 3-D ground electric field modeling such as the one applied by Marshalko et al. (2023) will be significantly less accurate and will not adequately represent the electric field in southern and central Sweden. Thus, we adopted the geoelectric field dataset described in Lanabere et al. (2023), which consists of 19 years (1 January 2000 to 30 December 2018) of 10-s geoelectric field data at 18 locations in Sweden covering six different latitudes and three longitudinal locations. The method used to obtain the dataset is described below.

2.1 IMAGE magnetic stations

In the previous study Lanabere et al. (2023) applied the 2-D SECS method (Vanhamäki and Juusola, 2020) with 10-s geomagnetic data from the IMAGE network to separate the measured B into the internal source (B_int) produced by induced telluric currents in the conducting ground, and the external source (B_ext) produced by electric currents in the ionosphere and magnetosphere. Due to the deficient spatial distribution of the IMAGE sites, the estimated B from the SECS method will be less accurate in Sweden than in Finland. Therefore, the authors used six IMAGE stations located along a similar longitude near Sweden (Glat 25.43°E), as shown with blue dots in Figure 1, and shifted B_ext to other locations (black dots), preferably to the same magnetic latitude. Geomagnetic data from the three Swedish IMAGE stations (pink dots) are used to estimate the method’s error.

Figure 1 (a) Location of the used IMAGE magnetic stations (blue and pink dots) and the corresponding shifted locations in Sweden (black dots) at the same geographic latitude, and the west, central, and east longitudes of the Swedish region at that latitude. Equal quasi-dipole geomagnetic latitudes (at 0 km) for the years 2000 and 2018 are presented along each station with blue and light-blue colors, respectively. (b–d) Horizontal distribution of the conductance S(x,y) for three sheets at different depths of thickness (D) developed by Korja et al. (2002). (b) The first layer (D1 = 0–10 km) corresponds to the upper crust, which includes the surface, seawater, and sediment. (c) The second layer (D2 = 10–30 km) corresponds to the middle crust. (d) The third layer (D3 = 30–60 km) corresponds to the lower crust. In all panels, the black dots indicate three equidistant locations (west, center, east) in Sweden at the same geographic latitude as the selected stations in color blue. Figure adapted from Lanabere et al. (2023).

2.2 Conductivity model

At each location shown in Figure 1a, the conductance is extracted as a function of depth from a three-layer conductance map for the Fennoscandian Shield (SMAP) developed by Korja et al. (2002), detailed further in Rosenqvist and Hall (2019). The conductance is divided into three layers of fixed depth and thickness: 1) the upper layer (0–10 km), 2) the middle layer (10–30 km), and 3) the lower layer (30–60 km). The locations and ground conductance (S) for each layer are shown in Figure 1. The largest conductance and longitudinal gradients are observed around the Lycksele magnetic station in the two upper layers. Lower conductance levels are observed in southwestern Sweden.

The geomagnetic field data is combined with ground electrical conductance into a 1-D conductivity model (Viljanen et al., 2004),

$$E_x\left(\omega \right)=\frac{Z\left(\omega \right)}{{\mu }_0}{B}_y\left(\omega \right), E_y\left(\omega \right)=-\frac{Z\left(\omega \right)}{{\mu }_0}{B}_x\left(\omega \right),$$(1)

where Z(ω) is the surface impedance of the layered Earth, defined by the thicknesses and electromagnetic parameters of the layers, which in this work is computed from the SMAP ground conductivity, μ₀ is the permeability of vacuum, and B_x and B_y are the north and east components of the local geomagnetic field, respectively. This approach is applied to 19 years of B_ext data at the various locations within Sweden to create an extensive dataset of geoelectric field data. This dataset is then utilized for the subsequent analysis.

3 Extreme value analysis

In this work, we apply the extreme value analysis to the 10-s geoelectric field data in Sweden. Rather than blocking the data into years and extracting the maximum from each block, it is more efficient to use the peak over threshold (POT) method (Coles, 2001). This method considers the raw data x = x₁,…,x_n a sequence of independent and identically distributed measurements and u a high threshold of x. Then, extreme events are identified by x_i:x_i > u. We label these exceedances by $x_{\left(1\right)},...,x_{\left(n_u\right)}$ which consist of n_u observations that exceed u, and define threshold excesses by y_j = x_(j)−u, for j = 1,...,n_u. The $\mathbf{y}=y_1,...,y_{n_u}$ may be regarded as independent realizations of a random variable whose distribution can be approximated by a member of the Generalized Pareto (GP) family. Then, the cumulative distribution function (G) of y is defined by

$$G_{\xi ,\sigma }\left(\mathbf{y}\right)=\{\begin{array}{ll}1-{\left(1+\frac{\xi \mathbf{y}}{\sigma }\right)}^{-\frac{1}{\xi }},\hspace{1em}& \mathrm{if}\hspace{1em}\mathrm{\xi }\ne 0\\ 1-\exp \left(-\frac{\mathbf{y}}{\sigma }\right),\hspace{1em}& \mathrm{if}\hspace{1em}\mathrm{\xi }=0,\end{array}$$(2)

where σ is the scale parameter, and ξ is the shape parameter. The threshold u is defined as the location parameter of the distribution. The sign of ξ determines the tail behavior. If ξ < 0, the distribution of excesses has an upper bound; if ξ = 0, the distribution is exponential, so unbounded; if ξ > 0 the distribution is heavy-tailed with no upper limit.

Within this study, x corresponds to the 10-s magnitude of the geoelectric field, and the threshold u corresponds to a high percentile of x. After a first inspection of the dataset, we observed that the data presents both time-dependent variation and short-term clustering. To illustrate this behavior, the exceedances over the 99th percentile (P99) of the |E| at the geographic latitude of station Muonio and the central longitude of Sweden are shown in Figure 2. The frequency of exceedances shows a temporal variation with the solar cycle, with more exceedances occurring during the declining phase of the solar cycle and fewer during solar minimum (2007–2010) (Fig. 2a). Moreover, it is evident that consecutive extreme values of |E| belong to the same cluster (Fig. 2b). Both consecutive exceedances and temporal variability with the solar cycle were expected since large perturbations of the ground magnetic field occur mainly during geomagnetic storms and severe geomagnetic storms are strongly modulated by the solar cycle (Chapman et al., 2020). Thus, it is crucial to properly deal with extreme event data points to obtain an independent exceedance dataset.

Figure 2 Example of |E| exceedance data showing time-dependent variation and short-term clustering. (a) Annual relative probability (Ci/N) of |E| exceedances where Ci is the number of elements in each year’s bin and N is the total number of elements (blue bars), compared to the monthly mean total sunspot number (red line) between 2000 and 2018. (b) Time series of exceedances from 11 October 2016 to 07 November 2016; the horizontal blue dashed line corresponds to the threshold. Both plots correspond to |E| data at the latitude of station Muonio (68.02°N) and central longitude in Sweden (20.52°E) at the station latitude, which exceeds the 99th percentile threshold of 0.12 V/km.

3.1 Declustering

If the series were independent, the threshold exceedances would be uniformly distributed. This is not observed in |E| (Figure 2). Thus, it is necessary to perform a declustering, which corresponds to a filtering of the dependent observations to obtain a set of threshold excesses that are approximately independent. In this study, we employed two different methods to define clusters of exceedances: the first is based on extremes separated by a fixed period of time, while the second method considers magnetospheric conditions to identify independent extremes. Then, we can assume that the cluster maxima are independent with conditional exceedances given by the GP (Eq. (2)).

In the first declustering method (DCM-A), we identified clusters by searching for extreme values that were separated by a fixed period of time. Consequently, the number of clusters depends on the selected threshold as well as the selected period of time. This declustering method has been applied before by several authors; for instance, Thomson et al. (2011) and Rogers et al. (2020) considered extreme values of 1-min geomagnetic measurements that exceed the 99.97% percentile to be part of the same cluster if the exceedances were separated by less than 12 h, while Meredith et al. (2016) considered a cluster to be active until the maxima dropped below the chosen threshold for a 9-h window. In this work, we establish that clusters should be separated by a period of at least three days without any exceedance. This approach of three days with low levels of |E| is intended to separate exceedances associated with different external drivers such as interplanetary coronal mass ejections. During solar cycles 23 and 24, their duration ranged from 3 to 90 h with an average of 28.8 h (Mishra et al., 2021). Upon reviewing the previous literature, we found a study by Hutchinson et al. (2011) where the authors analyzed 143 geomagnetic storms with SYM-H index and reported a mean duration of 2.3 days for weak storms, 2.8 days for moderate storms, and 3.9 days for intense storms. In these cases, at least 69% of the duration was part of the recovery phase, and smaller storms recovered more quickly than larger storms. Moreover, the larger values of |E| were observed near the main phase of the geomagnetic storm, where the largest magnetic field fluctuations occur. Thus, with a three-day separation between clusters, we intend to obtain independent observations, as it is longer than the aforementioned timescales. After identifying the clusters, only the maximum exceedance from each cluster is retained. An example of the DCM-A applied to |E| at the latitude of the Muonio station and central longitude, which exceeds the percentile threshold 99.0 and 99.7 is shown in Figure 3. In this period, we identify two clusters (gray shaded area) for exceedances over the 99.0 percentile (horizontal blue line), and the selected two maximum points for each cluster (blue dots). Similarly, it shows three clusters (dark gray area) of exceedances over the 99.7 percentile (horizontal orange line), and the three exceedance points for each cluster (blue and orange dots).

Figure 3 Example of declustering method A (DCM-A) being applied to the time series of |E| for the Muonio station’s latitude between 11 October 2016 and 07 November 2016. The gray shaded areas indicate clusters of exceedances over the 99.0th percentile (P99), while the dark gray shaded areas indicate clusters of exceedances over the 99.7th percentile (P99.7). Both are separated by at least three days. The horizontal blue and orange dashed lines correspond to the thresholds u = P99 and u = P99.7, respectively. The blue (blue and orange) points are the maximum exceedance inside each gray (dark gray) cluster.

A second declustering method (DCM-B) is employed since the first method does not ensure that the magnetosphere has returned to a quiet state. This may have two effects: a) we do not decluster if two consecutive short geomagnetic storms occur within 3 days, and b) we decluster a unique geomagnetic storm that is longer than 3 days. Therefore, we employ a relevant physical parameter to determine the activity state of the magnetosphere. First, we define the calm state of the magnetosphere based on the geomagnetic index SYM-H and the auroral index AE. This allows us to consider processes that drive geomagnetic storms and substorms. We used the lower threshold −50 nT < SYM-H as defined for a moderate geomagnetic storm by Gonzalez et al. (1994), and used the upper threshold of SYM-H < 15 nT as a limit for average SC/SI amplitude for magnetic clouds (29.01 nT) and noncloud ejectas (22.34 nT) (Veenadhari et al., 2012). Meanwhile, the calm substorm conditions are AE < 50 nT based on the previous work of Partamies et al. (2013). Second, we consider that two consecutive perturbed magnetospheric conditions are independent if they are separated by at least 6 h of calm magnetospheric conditions. This value is based on a previous study by Chu et al. (2015), in which the authors calculated the time between substorms, measuring the characteristic time of the unloading process in the magnetosphere. They found that the most probable waiting time between 1982 and 2012 is 80 min. The median and mean values are 3 and 8 h, respectively. Thus we choose a value twice the median value; the mean value is not considered because it is influenced by outliers.

As an example of the DCM-B, the defined calm (perturbed) magnetosphere conditions can be visualized in Figure 4a as white (gray) regions. The threshold limits for SYM-H and AE are shown as horizontal dashed lines. The number of clusters defined by this method does not depend on the selected threshold for E. There are a total of 1479 periods of perturbed magnetospheric conditions, and the longest cluster has a duration of 250 days from 27-Jan-2003 to 04-Oct-2003. This demonstrates why using a fixed time is not optimal since the active times with E not exceeding a high threshold can be much longer than three days. Then the maximum |E| exceedance from each cluster is retained; therefore, the number of clusters that contain an exceedance (i.e., the exceedance dataset) depends on the selected threshold. Figure 4b shows the exceedance with blue dots; in this example, the exceedance is the same for the P99 and P99.7 threshold selections.

Figure 4 Example of the declustering method B (DCM-B) applied to the time series of |E| for Muonio station latitude and central longitude in Sweden between 11 October 2016 and 07 November 2016.(a) Example of the time series of SYM-H and AE geomagnetic indices. The shaded gray area indicates perturbed magnetosphere periods that are separated by at least 6 h of calm magnetospheric conditions as defined by −50 nT < SYM-H < 15 nT and AE < 50 nT and are represented by white areas. (b) Time series of |E| from 11 October 2016 to 07 November 2016, the horizontal blue dashed line corresponds to the threshold value P99 and P99.7, and the blue dots correspond to the exceedance value in each cluster (gray shaded area).

The results of applying the declustering methods are compared with the same example of the unclustered distribution shown in Figure 2a. The solar cycle dependency of the exceedance data after applying DCM-A (Fig. 5a), shows a more uniform distribution compared to the unclustered distribution. However, for thresholds higher than 99.5 and at latitudes below 65°N, the distribution of exceedances still follows the solar cycle, with a low frequency of exceedances around the solar minimum between 2008 and 2011 (not shown). The distribution of the year of exceedance after applying DCM-B (shown in Fig. 5b) is more uniform compared to that in Figure 2a. The low frequency in 2003 resulted from the large cluster lasting 250 days. In conclusion, the exceedance dataset from which the Generalized Pareto distribution will be fitted and from which the return values will be computed depends on the declustering method and the selected threshold. Therefore, it is very important to choose the correct threshold. In this context, we applied both declustering methods to all six latitudes in Sweden and the three longitude locations (west, center, east), and for different thresholds between the 80th and 99.9th percentiles, with a step size of 1 between 80 and 99 and a step size of 0.1 between 99.0 and 99.9. Then, we determine the most appropriate threshold using two different methods, as described in the following subsection.

Figure 5 Example of the year of the |E| exceedance distribution (gray bars) and the monthly mean total sunspot number (red) between 2000 and 2018. (a) Exceedances after applying DCM-A, (b) Exceedances after applying DCM-B. Both plots correspond to |E| at the latitude of station Muonio (68.02°N) and the central longitude in Sweden (20.52°E) at the station latitude that exceeds the 99th percentile threshold.

3.2 Threshold selection

The optimal choice of the threshold should represent a balance between bias and variance. Setting a threshold too low is likely to violate the asymptotic motivation for the Generalized Pareto distribution, resulting in bias; conversely, a threshold set too high will yield too few excesses for accurate model estimation, contributing to increased variance (Coles, 2001). We used two different threshold selection methods: a) the parameter stability determined by fitting models across various thresholds, and b) the mean residual life (MRL) plot which is computed only from exceedance data. We applied both methods for each declustered data set across a range of different thresholds from the 80th percentile to the 99.9th percentile.

The parameter stability method states that if a GP distribution appropriately models the excesses above a threshold u₀, it follows that excesses above u (with u > u₀) should also follow a GP distribution with identical shape parameters for both distributions. Since the scale parameter changes with the threshold, a more appropriate quantity to compare distributions is the modified scale parameter computed as σ^* = σ_u−ξu (Coles, 2001). Then, above a level u₀ the estimates of both parameters ($\widehat{\xi }$, ${\widehat{\sigma }}^{\mathrm{*}}$) at a range of thresholds should be approximately constant. The 95% confidence interval for the estimated parameters are calculated as $C{I}_{\widehat{\xi }}=\widehat{\xi }\pm 1.96\sqrt{\mathrm{Var}\left(\widehat{\xi }\right)}$ using the variance-covariance matrix for $(\widehat{\sigma },\widehat{\xi })$ and $C{I}_{{\widehat{\sigma }}^{\mathrm{*}}}={\widehat{\sigma }}^{\mathrm{*}}\pm 1.96\sqrt{\mathrm{Var}\left({\widehat{\sigma }}^{\mathrm{*}}\right)}$ where $\mathrm{Var}\left({\widehat{\sigma }}^{\mathrm{*}}\right)$ is defined as:

$$\mathrm{Var}\left({\widehat{\sigma }}^{\mathrm{*}}\right)\approx \left[\begin{array}{ll}\frac{\partial {\sigma }^{\mathrm{*}}}{\partial {\sigma }_u},& \frac{\partial {\sigma }^{\mathrm{*}}}{\partial \xi }\end{array}\right]\left[\begin{array}{ll}\mathrm{Var}\left(\widehat{\sigma }\right)& \mathrm{Cov}(\widehat{\sigma },\widehat{\xi })\\ \mathrm{Cov}(\widehat{\xi },\widehat{\sigma })& \mathrm{Var}\left(\widehat{\xi }\right)\end{array}\right]\left[\begin{array}{l}\frac{\partial {\sigma }^{\mathrm{*}}}{\partial {\sigma }_u}\\ \frac{\partial {\sigma }^{\mathrm{*}}}{\partial \xi }\end{array}\right].$$(3)

We analyzed the function’s behavior within a specific range from the 95th percentile (P95) to the 99.9th percentile (P99.9). In this range, we considered windows of three to six consecutive values and calculated the variance for each window. Next, we fit a curve (such as a linear regression) to each window and computed the slope of the curve. We identified the window with the minimum variance, which suggests that the function remains the most constant within that interval. By reporting the slope of the fitted curve for this window, we can interpret how stable the function is; a slope close to zero indicates minimal change in the function across this region. This approach allows us to pinpoint the part of the function that exhibits the least variation.

The results of the threshold selection method for the exceedances using the DCM-A at all six latitudes and the central longitude in Sweden are shown in Figure 6. First, the shape parameter estimations are shown in Figure 6a–f. The linear fit for the five consecutive points window with the lowest variance is shown in cyan and the slope is included in each panel. The intervals with the lowest variance for other window sizes are shown in red. All latitudes, except for 62.25°N, exhibited the lowest variance between P99 and P99.5 across three different window sizes. We point out that the shape parameter at 62.25°N presents a notably higher value between P95 and P99 compared to P99.5 and P99.7. The modified scale parameter (Fig. 6g–l) shows large variations from P80 till P97 and the window with approximately constant values is observed between P99 and P99.5 at all latitudes, except at 62.25°N where it ranges from approximately P99.5 to P99.8. We note that at the two lowest latitudes the shape parameter presents the largest slope’s magnitude.

Figure 6 Parameter stability with DCM-A. Left: Shape parameter, Right: Modified scale parameter. Cyan lines represent the fitted line with the lowest variance, based on a window size of five consecutive points, with the slope reported in each panel. Red lines indicate the window with the lowest variance for alternative window sizes (three, four, and six consecutive points.

A similar analysis was done with the second exceedance dataset obtained using the declustering method DCM-B (Fig. 7). Both the shape and modified scale parameters present almost constant values between P96 and P99.4 at the three higher latitudes, represented by the low variance in this range considering the different window sizes and the low slope. However, at the two lower latitudes the range with lowest variance in the shape parameter ranges between P96 and P99.6 for different window sizes, and for latitude 62.25°N both shape and modified scale parameters present constant values for u > P99.4.

Figure 7 Parameter stability with DCM-B. Left: Shape parameter. Right: Modified scale parameter. Cyan lines represent the fitted line with the lowest variance, based on a window size of five consecutive points, with the slope reported in each panel. Red lines indicate the window with the lowest variance for alternative window sizes (three, four, and six consecutive points).

Based on the parameter stability method, both datasets indicate that the threshold among all latitudes should be selected between P99 and P99.5. However, due to the variable behavior within this range at latitude 62.25°N the threshold should be selected above P99.4. Since we want to compare the return values among the different locations, we aim to use the same threshold for all latitudes. Therefore, we selected P99.5 as the most suitable threshold. Next, we applied the second threshold selection method based on the mean residual life.

The MRL method is based on the mean of GP distribution. It states that if the GP distribution is valid for excess over a threshold u₀, it should be equally valid for all thresholds u > u₀. Then the GP mean value which is the same as the mean excess, is a linear function of u for u > u₀ (Coles, 2001). The MRL plot shows the relation between u and the mean of excess computed as $\frac{1}{n_u}{\mathrm{\Sigma }}_{j=1}^k\left(y_j\right)$ where y_j = x_i−u are the threshold excesses which consist of n_u observations that exceed u. The linearity is evaluated by applying a linear fit and calculating the normalized root mean square error NRMSE = RMSE /(y_max−y_min), where a value of zero would indicate a perfect fit to the data.

Figure 8 shows an example of the MRL plot with 95% confidence intervals using exceedance from both declustering methods (DCM-A and DCM-B) at all six latitudes and central longitude. We included the linear fit between P99.5 and P99.9 (red line) and the computed NRMSE. The mean excess plots using DCM-A (Fig. 8a–f) exhibit large values from the beginning of the plot (P80) and curve till P99 (purple dashed line), beyond which it is approximately linear (NRMSE < 0.12), except at 64.52N where the MRL presents an almost constant behavior and the NRMSE is the largest among all latitudes. We also applied the linear fit starting from P99 and calculated the NRMSE and Pearson correlation coefficients (r_p). The results are quite similar for thresholds P99 and P99.5 all returning values r_p > 0.94 and NRMSE < 0.12. This suggests that one should adopt a minimum value of u₀ = P99. The MRL plots with DCM-B (Fig. 8g–l) present a similar result to the one presented in DCM-A. There is a clear linear dependence on u beyond P99 (r_p > 0.94 and NRMSE < 0.10) at all latitudes except at 64.52 N. Moreover, both subsets show a negative MRL slope at the two higher latitudes, an almost zero slope at 64.52 N, and a positive MRL slope among the three lower latitudes. This behavior suggests that the tail of the exceedance distributions is similar among the mentioned latitudes.

Figure 8 Mean residual life plot for the |E| that exceed the threshold given by different percentiles between 80th and 99.9th (vertical dashed lines). The panels correspond to different latitudes in the central location in Sweden. Panels a–f: |E| exceedance using declustering method DCM-A. Panels g-l: |E| exceedance using declustering method DCM-B. The red line indicates the linear fit between thresholds P99.5 and P99.9. The normalized root mean square error is included in each panel.

After exploring both threshold selection methods, we conclude that at the two higher latitudes a threshold above P99 should be selected. At 64.52 N the MRL method presents the worst performance (largest NRMSE), but the parameter stability indicates a constant behavior of the parameters for values above the P99. At the three lower latitudes the MRL presents good performance for both P99 and P99.5 thresholds. However, the shape parameter is larger at P99 than at P99.5 and more evident at 62.25 N.

4 Results

In the previous Section we explained the steps to extract independent extreme events from a 19-years dataset of 10-s |E| estimated at 18 different locations in Sweden. Two exceedance data sets were obtained using two different cluster identifications: a) extreme events that are temporally separated by less than 3 days are considered to belong to the same cluster (DCM-A), and b) perturbed magnetospheric conditions that are separated by less 6 h of calms conditions are considered to belong to the same cluster (DCM-B). Then, for each method, the maximum |E| that exceed the 99.5th percentile in each cluster is retained for the final exceedance dataset. Finally the GP distribution can be fitted to each dataset.

4.1 Fitting the GP tail distribution

The maximum likelihood estimates of the Generalized Pareto parameters, together with 95% confidence intervals were determined to the exceedance |E| data over P99.5 at all locations declustered with methods DCM-A and DCM-B. The location, scale, and shape parameter results are shown in Figure 9. First, we analyzed the location parameter (i.e., the 99.5th percentile threshold), which is the same for both declustering methods. The threshold increases with latitude reaching its highest value where the largest |E| is obtained, except at geographic latitude 64.52 N at the center and east longitudes where the P99.5 is lower due to higher ground conductivity. This implies that higher latitudes present larger probabilities of large |E|. Another feature is that the location parameter suggests larger values at the west, at all latitudes due to lower conductance as shown in Lanabere et al. (2023). The larger difference among west, center, and east is observed at 64.52 N where large conductance gradients are present (see Fig. 1b–d). A similar behavior is observed in the scale parameter, mainly because σ changes with u (Coles, 2001). With regard to the declustering methods, σ presents lower values with DCM-A than DCM-B, which are more evident in west Sweden and 64.52 N.

Figure 9 Generalized Pareto distribution parameters (location, shape and scale parameters) obtained by fitting different exceedance |E| > P99.5 data sets. a) Exceedance set using DCM-A (dashed line) and DCM-B (solid line) at six different latitudes in Sweden and three longitudes (east, west center).

Finally, the shape parameter presents positive values at low latitudes (58.26 –62.25 N) while at the three higher latitudes ξ is positive close to zero, with lower 95% confidence limit within negative values. This behavior is similar to that obtained by Wintoft et al. (2016), except for latitude around 62 N where the authors present ξ values around zero. The difference may be attributed to the magnetic station’s location, as the authors selected a magnetic station in Scotland. This finding reinforces the results presented by Fogg et al. (2023) regarding the significance of more localized predictions. We highlight the behavior in the three parameters at 64.52 N. While the location and scale parameter show a large difference between west and east/center, the shape parameter remains essentially the same. This means that both u and ξ change with longitude due to different |E| related to changes in the conductance, but the behavior of the exceedance events remains the same. This characteristic remains the same between both declustering methods.

4.2 Return values

We are mainly interested in the largest |E| that is expected to be observed in N years, this value is known as the return value and was calculated following Coles (2001)

$$x_m=u+\frac{\sigma }{\xi }\left[{\left(m{\zeta }_u\theta \right)}^{\xi }-1\right],$$(4)

where m = N*n_y (n_y the number of observations per year, N number of years), σ and ξ are the scale and shape parameters of the Pareto distribution, ζ_u the probability of individual observations exceeding the threshold u (estimated as ${\widehat{\zeta }}_u=\frac{n_u}{n}$), and θ is the extremal index (estimated as $\widehat{\theta }=\frac{n_c}{n_u}$) with n the total number of data points, n_u the number of points exceeding the threshold u, and n_c is the number of clusters. We computed equation (4) at each location in Sweden (six latitudes and three longitudes) for different return periods (5, 10, 50, and 100 years).

The results are shown in Figure 10a–c for the exceedance set using DCM-B. The black crosses correspond to the 10 largest |E| values from the dataset at each latitude. At the three lower latitudes, the shape parameter is positive between 0.4 and 0.6, indicating that the distribution of extreme events have a tail which goes to zero slower than that of the three higher latitudes (ξ~0). As a consequence, at low latitudes the return values for long waiting times (50 and 100 years) can reach larger |E| compared with high latitudes. For example, at high latitudes the expected maximum in 100 years is between 2.0 and 2.5 V/km at the three different longitudes in Sweden. In central Sweden, the expected maximum in 100 years is between 4 and 6 V/km at the center and east longitudes and between 6 and 9 V/km in the west.

Figure 10 (a–c) Return values at west, center and east Sweden using DCM-B. The 10 largest modeled |E| values at each latitude are shown with black crosses. (d–f) Estimated return value normalized by the maximum |E| value, and 95% confidence interval at west, center, and east Sweden using DCM-B. The return levels at each latitude are artificially slightly shifted for a better visualization.

In our study, we analyzed 19 years of data, thus it is expected that we capture at least one event with a |E| around the 1-in-10 years event. Moreover, due to the distinct behavior of the tail, the estimated |E| values at high latitudes reach values similar to the expected 1-in-100 years, while at low latitudes, higher magnitudes are expected. When comparing the return values at different longitudes (i.e., different conductivity values) western Sweden exhibits larger expected values than the central and eastern regions. Interestingly, at 64.52°N the largest 1-in-100 event in central and eastern Sweden has a lower magnitude compared to the 1-in-1 year event in west Sweden. Furthermore, at lower latitudes, the maximum modeled |E| value is well below the expected value in 100 years. However, at higher latitudes the maximum |E| slightly exceeds the expected 1-in-100-years value, but still falls well within the limits of the error bar. The threshold selection methods performed robustly at these high latitudes, and the results remain consistent even when using a different threshold (P99). Thus, the 95% confidence intervals provide further interpretation of the results, offering a measure of the uncertainty around the estimated parameters.

Then, we included the 95% confidence interval for x_m derived by the delta method (Coles, 2001, page 82]. The 1 in 10- 50-, and 100 years return level plots normalized by the maximum |E| with 95% confidence intervals are shown in Figure 10d–f. The extremely large confidence intervals highlight the difficulties in modeling extreme values for long waiting times. In general, at lower latitudes, the maximum expected |E| in 100-years at least duplicates the current maximum value of the |E| and at higher latitudes the maximum expected is approximately the same order of magnitude as the maximum value modeled during the Halloween geomagnetic storm. The largest expected |E| in 100 years, given by the upper level of the confidence interval, is observed in western Sweden at 60.5°N. Its value is around five times the maximum modeled |E| value.

5 Validation of approximations

The results and analysis of the extreme events in Sweden were performed using E estimated from magnetic stations in Finland and Estonia, and ground conductivity in Sweden. In the previous study of Lanabere et al. (2023), the compilation of the full database was validated by comparing the distribution function of the daily maximum |E|, the validation presented here analyzes the distribution of the extreme events, which have a direct influence on the estimated return values. Thus, we performed the same analysis with the available data from the three Swedish stations that presented a long-time series. That is, we used the 10-s magnetic field measurements from 2000 to 2018 of the three a long in Sweden (KIR, LYC, UPS) to compute E in Sweden (without shifting) and using both B_ext and B_tot. We used the same declustering method and threshold as the previous results. That is, clusters of quiet and active magnetospheric conditions are defined by the AE and SYM-H geomagnetic indices. Then, we kept the maximum |E| value that exceeded the P99.5 threshold in each active cluster. Even though the three stations offer limited spatial coverage in Sweden such stations are used here to validate the use of the non-Swedish stations since it was shown in Lanabere et al. (2023) that this method likely underestimates slightly |E| for strong events.

First, we studied the effect of “shifting” B_ext. For this study we used B_ext data from the following pair of stations MUO/KIR, OUJ/LYC, and NUR/UPS (see Fig. 1) and the same conductance given by (KIR, LYC, and UPS respectively) to compute E. The GP fitting results applied to the three pairs of stations are shown with green and blue lines in Figure 11. We see that the location and scale parameters present a similar behavior. However, the shape parameter presents larger differences, where bigger values are obtained for the shifted B_ext at the two lower latitudes. A second approximation used in the estimation of E is related to the use of B_ext instead of B_tot (blue and pink lines, respectively). In this case, at all three latitudes the fitted parameters using B_tot are bigger than using B_ext but present the same relation among latitudes. The larger values in the location and scale parameters are expected since the larger E is obtained when using B_tot. Finally, the fitting parameters obtained from shifting B_ext to east Sweden are shown in the orange line (same as Fig. 9). The E used has both different B_ext and different but similar conductance. We see that the fitted values present the same behavior (tendency and sign) of previously fitted parameters. Thus, the two major approximations of the E estimation (B shifting, and B_ext) do not affect the behavior of the exceedance probability distribution function but may change the magnitude of the expected values, especially in south and central Sweden.

Figure 11 Generalized Pareto distribution parameters (location, shape, and scale parameters) obtained by fitting different exceedance |E| > P99.5 data sets using DCM-B with different magnetic data (B and conductance (S) for the E estimation. Orange: Bext shifted, S east Sweden, green: Bext shifted, S Sweden, Blue: (Bext shifted, S Sweden, Pink: Btot Sweden, S Sweden.

The estimated |E| and the return values for the six shifted stations to east Sweden and the three Swedish stations are shown in left panel of Figure 12. When comparing the 10 largest events (black crosses) we see that |E| presents similar values between the two closest locations. The largest differences are observed for long waiting times (50- and 100-years) around 60°N. The maximum E expected values at the three Swedish stations using B_ext and B_tot are shown in the right panel of Figure 12. As expected, larger |E| are obtained when using B_tot instead of B_ext at all latitudes. The large differences between |E_tot| and |E_ext| at 59.9°N for longer waiting times is due to the large shape parameter shown in Figure 11, which indicates that at this latitude the distribution of extreme |E_tot| events is slightly modified resulting in a heavier tail. For example at 59.9°N the expected value in 100 years is E_tot = 17.8 V/km and E_ext = 2.7 V/km. Meanwhile at high latitudes (67.84°N) E_tot = 3.1 V/km and E_ext = 2.1 V/km. At 64.61°N we notice that the maximum expected value in 100 years using B_ext is comparable with the maximum expected value in 1 year using B_tot. Thus, our results underestimate the expected values of E, and the difference depends on the latitude and on the waiting time. For this reason this results should be considered as a lower threshold for the return value.

Figure 12 Left: Return values east Sweden using (Finnish and Estonia) stations shifted to Sweden, exceedances over P99.5 and declustered using DCM-B (same as Fig. 10), and return values for the three Swedish stations. All return values were computed using Bext. The black crosses indicate the 10 largest |E| at each location. Right: Return values for the three Swedish stations using exceedances over P99.5 and declustered using DCM-B. Upward-pointing triangle symbol indicates that E was computed with Btot. Downward-pointing triangle symbol uses Bext (same as left panel).

6 Discussion

We performed an EV study in Sweden using 10-s resolution, and a national conductivity model that provided return values for different spatial locations in Sweden. We have benchmarked the return values taking into account the latitudinal changes in dB/dt as well as changes in |E| due to large scale changes in conductivity. Therefore this work is an important step in improving Sweden’s resilience to GIC events.

The purpose of this work is to study the behavior of the extreme geoelectric field in Sweden and present the expected return values in 10-, 50-, and 100-years which is a critical benchmarking and risk analysis that has not been done before. While previous studies have been conducted in Europe with the same intention, it is known that regional space weather predictions are needed on a national level. Because of the significant conductivity variations throughout Sweden, we cannot rely on previous studies for E values or rely on dB/dt. Due to the lack of spatial coverage of magnetometers in Sweden that present a long time series, such a localized study is not possible. However, in a recent study of Lanabere et al. (2023), the authors presented an alternative method to estimate the E in Sweden using the external part of the measured geomagnetic data from a close magnetic station combined with ground electrical conductance in Sweden. Thus, this study present an analysis of the extreme values of the 10-s |E| estimated for the period 2000 to 2018 at 18 different locations in Sweden.

One of the main challenges of the extreme value analysis is to create a sample of exceedances that are independent and equally distributed (Coles, 2001). Extreme values of E are usually clustered in time around the development of a geomagnetic storm and/or substorms. A widely used method to decluster data consists on looking for extreme values (based on the peaks over a high threshold) that are separated by a fixed period of time without any physical interpretation. To circumvent this issue, we devised a novel approach using physically relevant information about the state of the magnetosphere based on both geomagnetic storm and substorm indices (SYM-H and AE). The clusters of active magnetosphere are separated with periods of quiet conditions to ensure that the magnetosphere returns to a steady state and observations are independent. Then we constructed two extreme data sets at 18 different locations in Sweden. Two different methods were also applied to each dataset to select the most adequate threshold throughout different percentiles. Finally, the 99.5th percentile was used for the construction of the exceedance to estimate the Generalized Pareto parameters.

The results indicate that the distribution of extreme events presents a different behavior in latitude between North and South/Central Sweden. At higher latitudes (64.52°–68.02°N) the shape parameter is positive around zero indicating that the distribution falls faster to zero than at lower latitudes (58.26°–62.25°N) where the shape parameter (around 0.5) and its 95% confidence interval are positive. Similar shape values were obtained for the |E| at Uppsala (60.00°N) by Wintoft et al. (2016). This implies that the expected maximum |E| in 50- and 100-years are close to the maximum modeled |E| value in 2000–2018 at high latitudes, while at lower latitudes higher values are expected in the future or have been recorded in the past. Such as during the 13–14 July 1982 storm, when large |E| values (4–5 V/km) were estimated to occur in east Sweden around 59.4°N (Wik et al., 2009). A novel contribution from our results is that the maximum expected value in 100 years is likely to be located at 60.50°N. This value corresponds to three times the maximum |E| value, and up to 5 times larger in west Sweden when considering the 95th confidence interval. The latitudinal location of the largest expected |E| value aligns with the results of Myllys et al. (2014), where they found a maximum once in 100 years value between Solund and Karmøy (Glat = 59.21°–61.08°N) situated at a magnetic latitude similar to the geographic latitude 60.50°N in Sweden. The longitudinal differences in the distribution of the extremes rely on the location parameter, indicating that larger |E| occur in west Sweden than in central or east Sweden. This difference is larger at 64.52°N where a large horizontal gradient of the ground electrical conductivity is more pronounced.

Then we computed the E using the total and external part of the measured geomagnetic field from the three geomagnetic stations located in Sweden and reproduced the previous analysis and obtained the GP parameters. When using B_ext all parameters at the three latitudes fall within the 95% confidence intervals of the nearest location (using the same ground electrical conductivity), except for the shape parameter at Lycksele (64.61°N) which is smaller than in Oulujärvi shifted (64.52°N). The results using B_tot show larger parameters compared to using B_ext. The location and scale parameter differences are more evident with increasing latitude, while the shape parameter differences are bigger with decreasing latitude. As a result, the return values for short periods of time (1 in 1- and 5-years) are similar to those obtained with the shifted magnetic stations, but differences become larger with increasing waiting times.

7 Future implications

This study highlights that |E| in southern and central Sweden exhibit different behaviors during extreme space weather events with respect to northern Sweden. The EV analysis using dB/dt misses the effects of conductivity, which have been shown to be important, especially in regions with strong spatial gradients, such as at 64.52°N. Therefore, future studies should employ E instead of dB/dt to obtain useful return values that provide valuable information for GICs. Recently, on January 1, 2022 the IMAGE network included three geomagnetic stations in southern Sweden. Additionally, the recent strong geomagnetic storm during May 10–12, 2024 (estimated minimum Dst = −412 nT) will enable the inclusion of these new magnetic measurements in the SECS method for computing the E using 3-D modeling, and perform and new EV analysis in the future.

The mentioned statistical studies of the |E| in Sweden can be used to perform different risk scenarios in the whole Swedish electric grid. An example of GIC modeling studies in Sweden is the work performed by Wik et al. (2008) where the authors computed GIC in the 400 kV power grid in southern Sweden and compared with observations for 24 GIC events between 1998 and 2000. More recently Rosenqvist and Hall (2019) and Rosenqvist et al. (2022) developed a 3D modeling of GIC in Sweden which was compared and validated at two different transmission lines in north and south Sweden. Thus, the combination of the return values presented in this work with the 3D modeling of GIC developed by Rosenqvist et al. (2022) applied to the Swedish power grid can give a more accurate description of the nodes that are prone to GIC. This study can also be combined with recreated horizontal magnetic field at the surface from geospace simulations of past extreme events such as the Carrington event (Blake et al., 2021) or the response to an idealized interplanetary coronal mass ejection (Welling et al., 2021).

The mentioned latitudinal difference of the |E| is related mainly with the dynamic of ground magnetic perturbations (dB/dt) which relates with magnetosphere-ionosphere coupling processes (Wei et al., 2021). The auroral electrojet moves southward with enhanced geomagnetic activity, and it is the leading edge of the auroral electrojet that creates the largest geomagnetic variations. High latitudes experience the aurora daily when energy input from the Sun is low to medium. During extreme events, such as on May 10th, the auroral electrojet moves very far southward, resulting in lower magnetic field variations (and E field) at high latitudes. Long space missions such as the three Swarm spacecraft, the four spacecraft Cluster, the MMS constellation, and THEMIS constellation have opened the opportunity to conduct studies using simultaneous measurements and to go further in understanding the processes that connect ground, ionospheric, and magnetospheric processes (Aryan et al., 2022; Dong et al., 2023).

8 Summary and conclusions

We performed the extreme value analysis in Sweden using 10-s |E| at 18 locations from 2000 to 2018. The different locations in Sweden attempt to capture the changes in the ground conductivity. We meticulously selected extreme events to meet Generalized Pareto conditions and implemented two declustering and threshold selection methods. A new declustering method based on the magnetospheric conditions is defined and presented in this study. Comparing results with three Swedish stations validated our approach, revealing:

1. At higher geographic latitudes (64.52°–68.02°N) the distribution of extremes events falls faster to zero than at lower latitudes (58.26°–62.25°N) in Sweden. What has not been shown in previous studies is that at 60.50°N the distribution of extremes presents the heaviest tail.

2. The expected maximum |E| in 50- and 100 years are close to the current maximum value in the period 2000–2018 at high latitudes, while at lower latitudes higher values are expected in the future. A remarkable addition is that at lower latitudes in Sweden the 1 in 100 years return value normalized by the magnitude of the largest event corresponds to 2–3 times the current maximum value.

3. A finding not shown in previous studies is that larger return values are expected in the west than in center and east Sweden due to lower conductance in the west. A novel contribution of this study is the prediction that the largest expected maximum in 100-years is expected to occur around 60.50°N. An estimated lower threshold of the |E| at different locations in Sweden is summarized in Table 1. Larger values around 60°N could be obtained when considering the total magnetic field, with |E| increasing from 2.7 V/km when using B_ext to 17.8 V/km when using B_tot.

Acknowledgments

V. L., A. P. D., L. R and A. J. were funded by the Swedish Research Council project (Grant 2021-06259). A. P. D. received support from the Swedish National Space Agency (Grant 2020-00111). L.J. and A.V. were supported by Grant 339329 from the Academy of Finland. VL and AD thanks Adnane Osmane (University of Helsinki, Department of Physics, Helsinki, Finland) for helpful discussions. We thank the institutes who maintain the IMAGE Magnetometer Array: Tromsø Geophysical Observatory of UiT the Arctic University of Norway (Norway), Finnish Meteorological Institute (Finland), Institute of Geophysics Polish Academy of Sciences (Poland), GFZ German Research Centre for Geosciences (Germany), Geological Survey of Sweden (Sweden), Swedish Institute of Space Physics (Sweden), Sodankylä Geophysical Observatory of the University of Oulu (Finland), Polar Geophysical Institute (Russia), DTU Technical University of Denmark (Denmark), and Science Institute of the University of Iceland (Iceland). The provisioning of data from AAL, GOT, HAS, NRA, VXJ, FKP, ROE, BFE, BOR, HOV, SCO, KUL, and NAQ is supported by the ESA contracts number 4000128139/19/D/CT as well as 4000138064/22/D/KS. We thank the two anonymous reviewers for their insightful comments and suggestions, that helped improve the manuscript. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

Figure 1 (a) Location of the used IMAGE magnetic stations (blue and pink dots) and the corresponding shifted locations in Sweden (black dots) at the same geographic latitude, and the west, central, and east longitudes of the Swedish region at that latitude. Equal quasi-dipole geomagnetic latitudes (at 0 km) for the years 2000 and 2018 are presented along each station with blue and light-blue colors, respectively. (b–d) Horizontal distribution of the conductance S(x,y) for three sheets at different depths of thickness (D) developed by Korja et al. (2002). (b) The first layer (D1 = 0–10 km) corresponds to the upper crust, which includes the surface, seawater, and sediment. (c) The second layer (D2 = 10–30 km) corresponds to the middle crust. (d) The third layer (D3 = 30–60 km) corresponds to the lower crust. In all panels, the black dots indicate three equidistant locations (west, center, east) in Sweden at the same geographic latitude as the selected stations in color blue. Figure adapted from Lanabere et al. (2023).

In the text

Figure 2 Example of |E| exceedance data showing time-dependent variation and short-term clustering. (a) Annual relative probability (Ci/N) of |E| exceedances where Ci is the number of elements in each year’s bin and N is the total number of elements (blue bars), compared to the monthly mean total sunspot number (red line) between 2000 and 2018. (b) Time series of exceedances from 11 October 2016 to 07 November 2016; the horizontal blue dashed line corresponds to the threshold. Both plots correspond to |E| data at the latitude of station Muonio (68.02°N) and central longitude in Sweden (20.52°E) at the station latitude, which exceeds the 99th percentile threshold of 0.12 V/km.

In the text

Figure 3 Example of declustering method A (DCM-A) being applied to the time series of |E| for the Muonio station’s latitude between 11 October 2016 and 07 November 2016. The gray shaded areas indicate clusters of exceedances over the 99.0th percentile (P99), while the dark gray shaded areas indicate clusters of exceedances over the 99.7th percentile (P99.7). Both are separated by at least three days. The horizontal blue and orange dashed lines correspond to the thresholds u = P99 and u = P99.7, respectively. The blue (blue and orange) points are the maximum exceedance inside each gray (dark gray) cluster.

In the text

Figure 4 Example of the declustering method B (DCM-B) applied to the time series of |E| for Muonio station latitude and central longitude in Sweden between 11 October 2016 and 07 November 2016.(a) Example of the time series of SYM-H and AE geomagnetic indices. The shaded gray area indicates perturbed magnetosphere periods that are separated by at least 6 h of calm magnetospheric conditions as defined by −50 nT < SYM-H < 15 nT and AE < 50 nT and are represented by white areas. (b) Time series of |E| from 11 October 2016 to 07 November 2016, the horizontal blue dashed line corresponds to the threshold value P99 and P99.7, and the blue dots correspond to the exceedance value in each cluster (gray shaded area).

In the text

Figure 5 Example of the year of the |E| exceedance distribution (gray bars) and the monthly mean total sunspot number (red) between 2000 and 2018. (a) Exceedances after applying DCM-A, (b) Exceedances after applying DCM-B. Both plots correspond to |E| at the latitude of station Muonio (68.02°N) and the central longitude in Sweden (20.52°E) at the station latitude that exceeds the 99th percentile threshold.

In the text

	Figure 6 Parameter stability with DCM-A. Left: Shape parameter, Right: Modified scale parameter. Cyan lines represent the fitted line with the lowest variance, based on a window size of five consecutive points, with the slope reported in each panel. Red lines indicate the window with the lowest variance for alternative window sizes (three, four, and six consecutive points.
In the text

	Figure 7 Parameter stability with DCM-B. Left: Shape parameter. Right: Modified scale parameter. Cyan lines represent the fitted line with the lowest variance, based on a window size of five consecutive points, with the slope reported in each panel. Red lines indicate the window with the lowest variance for alternative window sizes (three, four, and six consecutive points).
In the text

Figure 8 Mean residual life plot for the |E| that exceed the threshold given by different percentiles between 80th and 99.9th (vertical dashed lines). The panels correspond to different latitudes in the central location in Sweden. Panels a–f: |E| exceedance using declustering method DCM-A. Panels g-l: |E| exceedance using declustering method DCM-B. The red line indicates the linear fit between thresholds P99.5 and P99.9. The normalized root mean square error is included in each panel.

In the text

	Figure 9 Generalized Pareto distribution parameters (location, shape and scale parameters) obtained by fitting different exceedance |E| > P99.5 data sets. a) Exceedance set using DCM-A (dashed line) and DCM-B (solid line) at six different latitudes in Sweden and three longitudes (east, west center).
In the text

	Figure 10 (a–c) Return values at west, center and east Sweden using DCM-B. The 10 largest modeled |E| values at each latitude are shown with black crosses. (d–f) Estimated return value normalized by the maximum |E| value, and 95% confidence interval at west, center, and east Sweden using DCM-B. The return levels at each latitude are artificially slightly shifted for a better visualization.
In the text

	Figure 11 Generalized Pareto distribution parameters (location, shape, and scale parameters) obtained by fitting different exceedance |E| > P99.5 data sets using DCM-B with different magnetic data (B and conductance (S) for the E estimation. Orange: Bext shifted, S east Sweden, green: Bext shifted, S Sweden, Blue: (Bext shifted, S Sweden, Pink: Btot Sweden, S Sweden.
In the text

Figure 12 Left: Return values east Sweden using (Finnish and Estonia) stations shifted to Sweden, exceedances over P99.5 and declustered using DCM-B (same as Fig. 10), and return values for the three Swedish stations. All return values were computed using Bext. The black crosses indicate the 10 largest |E| at each location. Right: Return values for the three Swedish stations using exceedances over P99.5 and declustered using DCM-B. Upward-pointing triangle symbol indicates that E was computed with Btot. Downward-pointing triangle symbol uses Bext (same as left panel).

In the text