Issue 
J. Space Weather Space Clim.
Volume 14, 2024



Article Number  22  
Number of page(s)  16  
DOI  https://doi.org/10.1051/swsc/2024025  
Published online  21 August 2024 
Research Article
Characterizing the distribution of extreme geoelectric field events in Sweden
^{1}
Swedish Institute of Space Physics, Institutet för rymdfysik, Box 537, 75121 Uppsala, Sweden
^{2}
Swedish Defense Research Agency, FOI Totalförsvarets forskningsinstitut, 164 90 Stockholm, Sweden
^{3}
Finnish Meteorological Institute, Dynamicum Erik Palménin aukio 1, 00560 Helsinki, Finland
^{*} Corresponding author: vanina.lanabere@irfu.se
Received:
19
February
2024
Accepted:
15
July
2024
Historically, Sweden has reported several impacts on transformers and transmission lines related to geomagnetically induced currents (GICs) that develop during strong space weather events. GICs are driven by the geoelectric field (E), and their intensity depends on various factors, including the lithology conductivity and the rate of change of the Earth’s magnetic field. The purpose of this study is to perform an extreme value (EV) analysis of the E magnitude at six different latitudes in Sweden and to express the maximum E that might be observed in 10, 50, and 100 years. We analyzed 10s E data in Sweden, obtained from a 1D model. This model incorporates 10s geomagnetic measurements from the IMAGE network and the vertical Earth’s ground electrical conductivity in Sweden, extracted from a 3D conductance map for the Fennoscandian region. Extreme E events tend to occur in clusters around geomagnetic disturbances (substorms and geomagnetic storms). Therefore, we applied two different methods to decluster the data. After declustering, Generalized Pareto (GP) distributions were fitted to the remaining extreme events that exceeded the 99.5th percentile. The EV analysis indicates that the shape parameter of the GP distribution depends on latitude. This implies that at higher geographic latitudes (64.52–68.02°N) the distribution decreases faster toward zero than at lower latitudes (58.26–62.25°N). As a result the expected maximum E in 100 years in central Sweden ranges between 4.0 and 8.5 V/km, while at higher latitudes, it ranges between 2.0 and 2.5 V/km, similar to the modeled geoelectric field values during the Halloween event in October 2003. In particular, around 60.50°N the distribution of extreme events exhibits the heaviest tail. When we additionally consider the effect of conductivity, the region of west Sweden around 60.50°N exhibits the largest expected maximum in 100 years with a value around 8.5 V/km. This is three times larger than the maximum modeled E at that latitude.
Key words: Geoelectric field / Extreme value theory / GIC
© V. Lanabere et al., Published by EDP Sciences 2024
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Extreme space weather events that represent a real threat to ground and spacebased infrastructure are categorized as lowfrequency but highconsequence 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/groundbased), physical properties of the hardware, time of the events, and the nature in which the external driver couples to nearEarth space are all important. For instance, during a strong solar event on 3–5 April 2010 the Intelsat satellite G15 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 highfrequency radio blackout that affected commercial airlines; excessive radiation exposure for highaltitude 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 groundbased 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 highvoltage 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, MaloneLeigh et al. (2024) generated geoelectric field hazard maps across Ireland using magnetotelluric transfer functions to model the geoelectric field over a 28year period using groundbased 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 stormtime 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 1min ground geomagnetic field data from 28 stations in various locations across Europe, including just two stations in Sweden. The authors employed the peakoverthreshold 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, 1min geomagnetic data from 14 stations, including two in Sweden, were analyzed using EV analysis applied to dB/dt and the 1D estimated E. The authors selected extreme events by taking the maximum values from oneyear 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 100years 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 highamplitude 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 lowresolution 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 1min 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 1s and 60s 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 10s 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 10s data and a highresolution 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 highresolution 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 100years 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 100years. 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 longtime 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 2D Spherical Elementary Current System (SECS) method (Vanhamäki and Juusola, 2020) which is used for 3D 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 10s 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 2D SECS method (Vanhamäki and Juusola, 2020) with 10s 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 quasidipole geomagnetic latitudes (at 0 km) for the years 2000 and 2018 are presented along each station with blue and lightblue 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 (D_{1} = 0–10 km) corresponds to the upper crust, which includes the surface, seawater, and sediment. (c) The second layer (D_{2} = 10–30 km) corresponds to the middle crust. (d) The third layer (D_{3} = 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 threelayer 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 1D conductivity model (Viljanen et al., 2004),
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, μ_{0} 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 10s 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_{1},…,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 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 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
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 heavytailed with no upper limit.
Within this study, x corresponds to the 10s 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 timedependent variation and shortterm 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 timedependent variation and shortterm clustering. (a) Annual relative probability (C_{i}/N) of E exceedances where C_{i} 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 (DCMA), 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 1min 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 9h 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 SYMH 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 threeday 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 DCMA 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 (DCMA) 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 (DCMB) 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 SYMH and the auroral index AE. This allows us to consider processes that drive geomagnetic storms and substorms. We used the lower threshold −50 nT < SYMH as defined for a moderate geomagnetic storm by Gonzalez et al. (1994), and used the upper threshold of SYMH < 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 DCMB, the defined calm (perturbed) magnetosphere conditions can be visualized in Figure 4a as white (gray) regions. The threshold limits for SYMH 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 27Jan2003 to 04Oct2003. 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 (DCMB) 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 SYMH 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 < SYMH < 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 DCMA (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 DCMB (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 DCMA, (b) Exceedances after applying DCMB. 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_{0}, it follows that excesses above u (with u > u_{0}) 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_{0} the estimates of both parameters (, ) at a range of thresholds should be approximately constant. The 95% confidence interval for the estimated parameters are calculated as using the variancecovariance matrix for and where is defined as:
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 DCMA 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 DCMA. 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 DCMB (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 DCMB. 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_{0}, it should be equally valid for all thresholds u > u_{0}. Then the GP mean value which is the same as the mean excess, is a linear function of u for u > u_{0} (Coles, 2001). The MRL plot shows the relation between u and the mean of excess computed as 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 (DCMA and DCMB) 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 DCMA (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_{0} = P99. The MRL plots with DCMB (Fig. 8g–l) present a similar result to the one presented in DCMA. 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 DCMA. Panels gl: E exceedance using declustering method DCMB. 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 19years dataset of 10s 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 (DCMA), and b) perturbed magnetospheric conditions that are separated by less 6 h of calms conditions are considered to belong to the same cluster (DCMB). 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 DCMA and DCMB. 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 DCMA than DCMB, 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 DCMA (dashed line) and DCMB (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)
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 ), and θ is the extremal index (estimated as ) 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 DCMB. 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 DCMB. 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 DCMB. 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 1in10 years event. Moreover, due to the distinct behavior of the tail, the estimated E values at high latitudes reach values similar to the expected 1in100 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 1in100 event in central and eastern Sweden has a lower magnitude compared to the 1in1 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 1in100years 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 100years 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 longtime series. That is, we used the 10s 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 SYMH 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 nonSwedish 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 DCMB with different magnetic data (B and conductance (S) for the E estimation. Orange: B_{ext} shifted, S east Sweden, green: B_{ext} shifted, S Sweden, Blue: (B_{ext} shifted, S Sweden, Pink: B_{tot} 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 100years) 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 DCMB (same as Fig. 10), and return values for the three Swedish stations. All return values were computed using B_{ext}. 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 DCMB. Upwardpointing triangle symbol indicates that E was computed with B_{tot}. Downwardpointing triangle symbol uses B_{ext} (same as left panel). 
6 Discussion
We performed an EV study in Sweden using 10s 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 100years 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 10s 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 (SYMH 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 100years 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 5years) 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 3D 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 magnetosphereionosphere 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 10s 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:
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.
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.
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 100years 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}.
Expected maximum E in V/km at different locations in Sweden for different waiting times.
Acknowledgments
V. L., A. P. D., L. R and A. J. were funded by the Swedish Research Council project (Grant 202106259). A. P. D. received support from the Swedish National Space Agency (Grant 202000111). 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
 Allen J. 2010. The Galaxy 15 anomaly: another satellite in the wrong place at a critical time. Space Weather 8(6): S06008. https://doi.org/10.1029/2010SW000588. [Google Scholar]
 Alves Ribeiro J, Pinheiro FJG, Pais MA. 2021. First estimations of geomagnetically induced currents in the South of Portugal. Space Weather 19(1): e2020SW002, 546. https://doi.org/10.1029/2020SW002546. [CrossRef] [Google Scholar]
 Alves Ribeiro J, Pinheiro FJG, Pais MA, Santos R, Cardoso J, BaltazarSoares P, Monteiro Santos FA. 2023. Toward more accurate GIC estimations in the Portuguese power network. Space Weather 21(6): e2022SW003, 397. https://doi.org/10.1029/2022SW003397. [CrossRef] [Google Scholar]
 Aryan H, Bortnik J, Li J, Weygand JM, Chu X, Angelopoulos V. 2022. Multiple conjugate observations of magnetospheric fast flow bursts using THEMIS observations. Ann Geophys 40(4): 531–544. https://doi.org/10.5194/angeo405312022. [CrossRef] [Google Scholar]
 Bergin A, Chapman SC, Watkins NW, Moloney NR, Gjerloev JW. 2023. Extreme event statistics in Dst, SYMH, and SMR geomagnetic indices. Space Weather 21(3): e2022SW003, 304. https://doi.org/10.1029/2022SW003304. [CrossRef] [Google Scholar]
 Blake SP, Pulkkinen A, Schuck PW, Glocer A, Oliveira DM, Welling DT, Weigel RS, Quaresima G. 2021. Recreating the horizontal magnetic field at Colaba during the Carrington event with geospace simulations. Space Weather 19(5): e2020SW002, 585. https://doi.org/10.1029/2020SW002585. [CrossRef] [Google Scholar]
 Chapman SC, McIntosh SW, Leamon RJ, Watkins NW. 2020. Quantifying the solar cycle modulation of extreme space weather. Geophys Res Lett 47(11): e2020GL087, 795. https://doi.org/10.1029/2020GL087795. [CrossRef] [Google Scholar]
 Chu X, McPherron RL, Hsu TS, Angelopoulos V. 2015. Solar cycle dependence of substorm occurrence and duration: implications for onset. J Geophys Res Space Phys 120(4): 2808–2818. https://doi.org/10.1002/2015JA021104. [CrossRef] [Google Scholar]
 Cliver EW, Dietrich WF. 2013. The 1859 space weather event revisited: limits of extreme activity. J Space Weather Space Clim 3: A31. https://doi.org/10.1051/swsc/2013053. [Google Scholar]
 Coles S. 2001. An introduction to statistical modeling of extreme values. Springer Series in Statistics, Springer, London. ISBN 9781447136750. https://doi.org/10.1007/9781447136750. [CrossRef] [Google Scholar]
 Dimmock AP, Rosenqvist L, Welling DT, Viljanen A, Honkonen I, Boynton RJ, Yordanova E. 2020. On the regional variability of dB/dt and its significance to GIC. Space Weather 18(8): e2020SW002, 497. https://doi.org/10.1029/2020SW002497. [CrossRef] [Google Scholar]
 Doherty P, Coster A, Murtagh W. 2004. Space weather effects of October–November 2003. GPS Sol 8(4): 267–271. https://doi.org/10.1007/s1029100401093. [CrossRef] [Google Scholar]
 Dong XC, Dunlop MW, Xiao C, Wei D, Wang TY, Zhao JS. 2023. Simultaneous mesoscale polar cusp fieldaligned currents measured on mid and lowaltitude satellites. Geophys Res Lett 50(1): e2022GL102, 460. https://doi.org/10.1029/2022GL102460. [CrossRef] [Google Scholar]
 Eastwood JP, Biffis E, Hapgood MA, Green L, Bisi MM, Bentley RD, Wicks R, McKinnell LA, Gibbs M, Burnett C. 2017. The economic impact of space weather: where do we stand? Risk Anal 37(2): 206–218. https://doi.org/10.1111/risa.12765. [CrossRef] [Google Scholar]
 Fang TW, Kubaryk A, Goldstein D, Li Z, FullerRowell T, Millward G, Singer HJ, Steenburgh R, Westerman S, Babcock E. 2022. Space weather environment during the SpaceX starlink satellite loss in February 2022. Space Weather 20(11): e2022SW003, 193. https://doi.org/10.1029/2022SW003193. [CrossRef] [Google Scholar]
 Fogg AR, Jackman CM, MaloneLeigh J, Gallagher PT, Smith AW, Lester M, Walach MT, Waters JE. 2023. Extreme value analysis of ground magnetometer observations at Valentia observatory, Ireland. Space Weather 21(7): e2023SW003, 565. https://doi.org/10.1029/2023SW003565. [CrossRef] [Google Scholar]
 Gannon JL, Birchfield AB, Shetye KS, Overbye TJ. 2017. A comparison of peak electric fields and GICs in the Pacific Northwest using 1D and 3D conductivity. Space Weather 15(11): 1535–1547. https://doi.org/10.1002/2017SW001677. [CrossRef] [Google Scholar]
 Gaunt CT, Coetzee G. 2007. Transformer failures in regions incorrectly considered to have low GICrisk. In: 2007 IEEE Lausanne Power Tech, pp. 807–812. https://doi.org/10.1109/PCT.2007.4538419. [CrossRef] [Google Scholar]
 Gonzalez WD, Joselyn JA, Kamide Y, Kroehl HW, Rostoker G, Tsurutani BT, Vasyliunas VM. 1994. What is a geomagnetic storm? J Geophys Res Space Phys 99(A4): 5771–5792. https://doi.org/10.1029/93JA02867. [CrossRef] [Google Scholar]
 Hapgood M. 2019. The great storm of May 1921: an exemplar of a dangerous space weather event. Space Weather 17(7): 950–975. https://doi.org/10.1029/2019SW002195. [NASA ADS] [CrossRef] [Google Scholar]
 Hutchinson JA, Wright DM, Milan SE. 2011. Geomagnetic storms over the last solar cycle: a superposed epoch analysis. J Geophys Res Space Phys 116(A9). https://doi.org/10.1029/2011JA016463. [Google Scholar]
 Kappenman JG. 2006. Great geomagnetic storms and extreme impulsive geomagnetic field disturbance events – an analysis of observational evidence including the great storm of May 1921. Adv Space Res 38(2): 188–199. The Great Historical Geomagnetic Storm of 1859: A Modern Look. https://doi.org/10.1016/j.asr.2005.08.055. [CrossRef] [Google Scholar]
 Kataoka R, Shiota D, Fujiwara H, Jin H, Tao C, Shinagawa H, Miyoshi Y. 2022. Unexpected space weather causing the reentry of 38 Starlink satellites in February 2022. J Space Weather Space Clim 12: 41. https://doi.org/10.1051/swsc/2022034. [CrossRef] [EDP Sciences] [Google Scholar]
 Koons HC. 2001. Statistical analysis of extreme values in space science. J Geophys Res Space Phys 106(A6): 10915–10921. https://doi.org/10.1029/2000JA000234. [CrossRef] [Google Scholar]
 Korja T, Engels M, Zhamaletdinov AA, Kovtun AA, Palshin NA, Smirnov MY, Tokarev AD, Asming VE, Vanyan LL, Vardaniants IL. 2002. Crustal conductivity in Fennoscandia—a compilation of a database on crustal conductance in the Fennoscandian Shield. Earth Planets Space 54: 535–558. https://doi.org/10.1186/BF03353044. [Google Scholar]
 Lanabere V, Dimmock AP, Rosenqvist L, Juusola L, Viljanen A, Johlander A, Odelstad E. 2023. Analysis of the geoelectric field in Sweden over solar cycles 23 and 24: spatial and temporal variability during strong GIC events. Space Weather 21(12): e2023SW003, 588. https://doi.org/10.1029/2023SW003588. [CrossRef] [Google Scholar]
 Loto’aniu TM, Singer HJ, Rodriguez JV, Green J, Denig W, Biesecker D, Angelopoulos V. 2015. Space weather conditions during the Galaxy 15 spacecraft anomaly. Space Weather 13(8): 484–502. https://doi.org/10.1002/2015SW001239. [CrossRef] [Google Scholar]
 Love JJ. 2021. Extremeevent magnetic storm probabilities derived from rank statistics of historical Dst intensities for solar cycles 14–24. Space Weather 19(4): e2020SW002, 579. https://doi.org/10.1029/2020SW002579. [CrossRef] [Google Scholar]
 Love JJ, Lucas GM, Rigler EJ, Murphy BS, Kelbert A, Bedrosian PA. 2022. Mapping a magnetic superstorm: March 1989 geoelectric hazards and impacts on United States power systems. Space Weather 20(5): e2021SW003, 030. https://doi.org/10.1029/2021SW003030. [Google Scholar]
 MaloneLeigh J, Campanyà J, Gallagher PT, Hodgson J, Hogg C. 2024. Mapping geoelectric field hazards in Ireland. Space Weather 22(2): e2023SW003, 638. https://doi.org/10.1029/2023SW003638. [CrossRef] [Google Scholar]
 Marshalko E, Kruglyakov M, Kuvshinov A, Viljanen A. 2023. Threedimensional modeling of the ground electric field in Fennoscandia during the Halloween geomagnetic storm. Space Weather 21(9): e2022SW003, 370. https://doi.org/10.1029/2022SW003370. [CrossRef] [Google Scholar]
 Meredith NP, Horne RB, Isles JD, Green JC. 2016. Extreme energetic electron fluxes in low Earth orbit: analysis of POES E > 30, E > 100, and E > 300 keV electrons. Space Weather 14(2): 136–150. https://doi.org/10.1002/2015SW001348. [CrossRef] [Google Scholar]
 Meredith NP, Horne RB, Isles JD, Rodriguez JV. 2015. Extreme relativistic electron fluxes at geosynchronous orbit: analysis of GOES E > 2 MeV electrons. Space Weather 13(3): 170–184. https://doi.org/10.1002/2014SW001143. [CrossRef] [Google Scholar]
 Meredith NP, Horne RB, Sandberg I, Papadimitriou C, Evans HDR. 2017. Extreme relativistic electron fluxes in the Earth’s outer radiation belt: analysis of INTEGRAL IREM data. Space Weather 15(7): 917–933. https://doi.org/10.1002/2017SW001651. [CrossRef] [Google Scholar]
 Mishra W, Doshi U, Srivastava N. 2021. Radial sizes and expansion behavior of ICMEs in solar cycles 23 and 24. Front Astron Space Sci 8: 142. https://doi.org/10.3389/fspas.2021.713999. [CrossRef] [Google Scholar]
 Myllys M, Viljanen A, Rui OA, Ohnstad TM. 2014. Geomagnetically induced currents in Norway: the northernmost highvoltage power grid in the world. J Space Weather Space Clim 4: A10. https://doi.org/10.1051/swsc/2014007. [CrossRef] [EDP Sciences] [Google Scholar]
 Partamies N, Juusola L, Tanskanen E, Kauristie K. 2013. Statistical properties of substorms during different storm and solar cycle phases. Ann Geophys 31(2): 349–358. https://doi.org/10.5194/angeo313492013. [CrossRef] [Google Scholar]
 Patterson CJ, Wild JA, Boteler DH. 2023. Modeling the impact of geomagnetically induced currents on electrified railway signaling systems in the United Kingdom. Space Weather 21(3): e2022SW003, 385. https://doi.org/10.1029/2022SW003385. [CrossRef] [Google Scholar]
 Pulkkinen A, Bernabeu E, Thomson A, Viljanen A, Pirjola R, et al. 2017. Geomagnetically induced currents: science, engineering, and applications readiness. Space Weather 15(7): 828–856. https://doi.org/10.1002/2016SW001501. [NASA ADS] [CrossRef] [Google Scholar]
 Pulkkinen A, Lindahl S, Viljanen A, Pirjola R. 2005. Geomagnetic storm of 29–31 October 2003: geomagnetically induced currents and their relation to problems in the Swedish highvoltage power transmission system. Space Weather 3(8). https://doi.org/10.1029/2004SW000123. [CrossRef] [Google Scholar]
 Pulkkinen A, Viljanen A, Pirjola R. 2006. Estimation of geomagnetically induced current levels from different input data. Space Weather 4(8). https://doi.org/10.1029/2006SW000229. [CrossRef] [Google Scholar]
 Rogers NC, Wild JA, Eastoe EF, Gjerloev JW, Thomson AWP. 2020. A global climatological model of extreme geomagnetic field fluctuations. J Space Weather Space Clim 10: 5. https://doi.org/10.1051/swsc/2020008. [CrossRef] [EDP Sciences] [Google Scholar]
 Rosenqvist L, Fristedt T, Dimmock AP, Davidsson P, Fridström R, et al. 2022. 3D modeling of geomagnetically induced currents in Sweden—validation and extreme event analysis. Space Weather 20(3): e2021SW002, 988. https://doi.org/10.1029/2021SW002988. [CrossRef] [Google Scholar]
 Rosenqvist L, Hall JO. 2019. Regional 3D modeling and verification of geomagnetically induced currents in Sweden. Space Weather 17(1): 27–36. https://doi.org/10.1029/2018SW002084. [CrossRef] [Google Scholar]
 Schillings A, Palin L, Opgenoorth HJ, Hamrin M, Rosenqvist L, Gjerloev JW, Juusola L, Barnes R. 2022. Distribution and occurrence frequency of dB/dt spikes during magnetic storms 1980–2020. Space Weather 20(5): e2021SW002, 953. https://doi.org/10.1029/2021SW002953. [CrossRef] [Google Scholar]
 Thomson AWP, Dawson EB, Reay SJ. 2011. Quantifying extreme behavior in geomagnetic activity. Space Weather 9(10). https://doi.org/10.1029/2011SW000696. [Google Scholar]
 Trivedi NB, Vitorello I, Kabata W, Dutra SLG, Padilha AL, et al. 2007. Geomagnetically induced currents in an electric power transmission system at low latitudes in Brazil: a case study. Space Weather 5(4). https://doi.org/10.1029/2006SW000282. [CrossRef] [Google Scholar]
 Tsubouchi K, Omura Y. 2007. Longterm occurrence probabilities of intense geomagnetic storm events. Space Weather 5(12). https://doi.org/10.1029/2007SW000329. [CrossRef] [Google Scholar]
 Vanhamäki H, Juusola L. 2020. Introduction to spherical elementary current systems, 5–33. Springer International Publishing, Cham. ISBN 9783030267322. https://doi.org/10.1007/9783030267322_2. [Google Scholar]
 Veenadhari B, Selvakumaran R, Singh R, Maurya AK, Gopalswamy N, Kumar S, Kikuchi T. 2012. Coronal mass ejection–driven shocks and the associated sudden commencements/sudden impulses. J Geophys Res Space Phys 117(A4). https://doi.org/10.1029/2011JA017216. [CrossRef] [Google Scholar]
 Viljanen A, Pulkkinen A, Amm O, Pirjola R, Korja T, BEAR Working Group. 2004. Fast computation of the geoelectric field using the method of elementary current systems and planar Earth models. Ann Geophys 22(1): 101–113. https://doi.org/10.5194/angeo221012004. [CrossRef] [Google Scholar]
 Watari S, Kunitake M, Kitamura K, Hori T, Kikuchi T, et al. 2009. Measurements of geomagnetically induced current in a power grid in Hokkaido, Japan. Space Weather 7(3). https://doi.org/10.1029/2008SW000417. [Google Scholar]
 Webb DF, Allen JH. 2004. Spacecraft and ground anomalies related to the October–November 2003 solar activity. Space Weather 2(3). https://doi.org/10.1029/2004SW000075. [Google Scholar]
 Wei D, Dunlop MW, Yang J, Dong X, Yu Y, Wang T. 2021. Intense dB/dt variations driven by nearEarth Bursty Bulk Flows (BBFs): a case study. Geophys Res Lett 48(4): e2020GL091, 781. https://doi.org/10.1029/2020GL091781. [CrossRef] [Google Scholar]
 Welling DT, Love JJ, Rigler EJ, Oliveira DM, Komar CM, Morley SK. 2021. Numerical simulations of the geospace response to the arrival of an idealized perfect interplanetary coronal mass ejection. Space Weather 19(2): e2020SW002, 489. https://doi.org/10.1029/2020SW002489. [CrossRef] [Google Scholar]
 Wik M, Pirjola R, Lundstedt H, Viljanen A, Wintoft P, Pulkkinen A. 2009. Space weather events in July 1982 and October 2003 and the effects of geomagnetically induced currents on Swedish technical systems. Ann Geophys 27(4): 1775–1787. https://doi.org/10.5194/angeo2717752009. [CrossRef] [Google Scholar]
 Wik M, Viljanen A, Pirjola R, Pulkkinen A, Wintoft P, Lundstedt H. 2008. Calculation of geomagnetically induced currents in the 400 kV power grid in southern Sweden. Space Weather 6(7): 07005. https://doi.org/10.1029/2007SW000343. [Google Scholar]
 Wintoft P, Viljanen A, Wik M. 2016. Extreme value analysis of the time derivative of the horizontal magnetic field and computed electric field. Ann Geophys 34(4): 485–491. https://doi.org/10.5194/angeo344852016. [CrossRef] [Google Scholar]
 Xue D, Yang J, Liu Z, Yu S. 2023. Examining the economic costs of the 2003 Halloween storm effects on the North Hemisphere aviation using flight data in 2019. Space Weather 21(3): e2022SW003, 381. https://doi.org/10.1029/2022SW003381. [CrossRef] [Google Scholar]
Cite this article as: Lanabere V, Dimmock AP, Rosenqvist L, Viljanen A, Juusola L, et al. 2024. Characterizing the distribution of extreme geoelectric field events in Sweden. J. Space Weather Space Clim. 14, 22. https://doi.org/10.1051/swsc/2024025.
All Tables
Expected maximum E in V/km at different locations in Sweden for different waiting times.
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 quasidipole geomagnetic latitudes (at 0 km) for the years 2000 and 2018 are presented along each station with blue and lightblue 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 (D_{1} = 0–10 km) corresponds to the upper crust, which includes the surface, seawater, and sediment. (c) The second layer (D_{2} = 10–30 km) corresponds to the middle crust. (d) The third layer (D_{3} = 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 timedependent variation and shortterm clustering. (a) Annual relative probability (C_{i}/N) of E exceedances where C_{i} 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 (DCMA) 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 (DCMB) 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 SYMH 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 < SYMH < 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 DCMA, (b) Exceedances after applying DCMB. 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 DCMA. 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 DCMB. 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 DCMA. Panels gl: E exceedance using declustering method DCMB. 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 DCMA (dashed line) and DCMB (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 DCMB. 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 DCMB. 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 DCMB with different magnetic data (B and conductance (S) for the E estimation. Orange: B_{ext} shifted, S east Sweden, green: B_{ext} shifted, S Sweden, Blue: (B_{ext} shifted, S Sweden, Pink: B_{tot} 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 DCMB (same as Fig. 10), and return values for the three Swedish stations. All return values were computed using B_{ext}. 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 DCMB. Upwardpointing triangle symbol indicates that E was computed with B_{tot}. Downwardpointing triangle symbol uses B_{ext} (same as left panel). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.