Interactive comment on “ Defining scale thresholds for geomagnetic storms through statistics

The paper attempts to define some threshold levels from statistics in several magnetic indices. The initial idea was attractive because it promised to give some values which could be taken as more “objective” respect the usual ones. However, the results fail in the sense that the resulting thresholds values are disconnected from the physical causes that origin those disturbances either from the consequences provoked by them, which, in fact, are the purpose of this journal. Moreover, defining the threshold as the intersection of the cumulative of a model distribution and the data distribution looks more a mathematical artifact than a fact. To take advantage of the work done, I’d suggest to the authors to submit this paper to another more academic journal changing

Quantifying geomagnetic storm intensities ::::::::::: disturbances : is one of the main aims in space weather.Extreme space weather events are treated now as natural environmental hazards since they pose a serious threat to different technologies, including critical infrastructures.The threats include power grid upsets, that can be damaged due to sudden changes in the Earth's magnetic field, originating :::::: leading :: to : geomagnetically induced currents (GICs); problems in GNSS positioning; orbit displacement of satellites, or radio blackouts, among other undesirable effects.A summary of these risks and their economic and societal impact appeared recently in Eastwood et al. (2017).Very recently, macro-economic losses have been modelled considering worst-case scenarios centered in different areas of the globe (Schulte in den Bäumen et al., 2014).
Different geomagnetic indices have been created to account for different geophysical phenomena, and then, a scale was defined accordingly.Scales may rely on the physics behind an index, or a scale may be effect-based only.Some examples of the kind of scales are for earthquakes -effect-based Mercalli scale, and the M w scale, that is physics-based; for nuclear risks, the INES scale is an effect-based, and the NAMS scale is more physical (Wheatley et al., 2017).Below, we comment the most commonly employed indices for space weather.
The Dst is a geomagnetic index, considered to provide a quantitative measurement of geomagnetic disturbances at low latitudes.The method for the elaboration of the Dst index is described in :: by :::::::::::::::::::::::: Sugiura and Kamei (1991) : , and it is assumed as a representation of the ring current enhancement encircling the Earth.This index is widely used in space weather.However, the Dst index has been questioned to misrepresent the ring current, not accounting for asymmetries (Grafe, 1998); or being inaccurate due to important seasonal offsets (Karinen and Mursula, 2006), that may even lead to a mismatch in the currently used scale.
All this variety may pose an additional complication for the user to differentiate between physical and effect-based scales, their thresholds and their applicability range.In the end, proper threshold definition is necessary for considerations on the importance of geomagnetic disturbances and for further utilization by users and decision makers.
In this paper, we use different distribution functions seeking the best fit for the regional index LDiñ and its derivative LCiñ; then, we define thresholds for this index and its derivative.Fitting distributions aiming at recurrence or waiting times estimations is out of the scope of this work.Opposite to the vast majority of the mentioned references, we do not use thresholds to help fitting a part or the whole distribution; threshold rationale and distribution properties are the main purposes.Hence, one of the key goals in this paper is aiming at overcoming some of the arbitrariness to establish thresholds for geomagnetic disturbances and its rate of change through the objective method explained in this paper.Its application can be extended to any other global or regional index.
The paper is structured as follows: data used are presented in Section §2, along with the preliminary statistical considerations; next in Section §3, the method is described, obtaining the best fit distributions to the indices in Section §3.1; and then in Section §3.2, the thresholds for these indices are defined by these statistical best distribution intersects.The results and discussion are presented in Section §4, and Conclusions appear in Section §5.

Geomagnetic data
We used LDiñ data computed from SPT geomagnetic observatory (San Pablo de los Montes, Toledo, Spain) from 1997 to 2012 -solar cycles 23 and 24, comprising two rising phases and one solar maximum -, processed in such a way to remove the daily variation in the H component of the geomagnetic field, which is not straightforward for mid-latitude data.The processing method, which is now patent pending (Guerrero et al., 2016), will be actually detailed in a forthcoming paper.
Therefore we consider LDiñ as this local geomagnetic index after this specific processing, whose units are nT and its temporal cadence is 1 minute.We define hereafter LCiñ as the temporal derivative of LDiñ, in nT min −1 .Since the goal of this paper is not the storm recurrence or waiting time, not timing between peak values is considered.The whole dataset are 1-min values without any further consideration of recurrence. :: To :::::: prove ::: the :::::::: validity :: of ::::: LDiñ ::: for :::::::::: monitoring ::::: local :::::::::::: geomagnetic ::::::::::: disturbances :::::: while :::::::::: effectively ::::::::: removing ::::: other ::::::: trends, ::: we :::::  To give a hint about the value distributions of LDiñ and LCiñ, we plotted the histograms, normalized to the total amount, in Fig. 3.The bin size is 5 nT for LDiñ and 0.5 nT min −1 bin size for LCiñ.On geomagnetic storms, negative values actually define the storm since they mean a horizontal component depression of the geomagnetic field due to a ring current enhancement; therefore we will take into account only these negative values, naming it as |neg(LDiñ)|.It is evident that this distribution is not gaussian but other kind of distribution, with a long tail and very skewed to negative values, as shown in in Fig. 3 (top left).
We consider that the positive and negative values of LCiñ are meaningful because both may be related to the geomagnetically induced current (GIC) generation, interesting for space weather purposes.Therefore, we will consider this possibility equally for both signs by studying the absolute value of LCiñ, naming it |LCiñ|.As shown in Fig. 3  3 Distribution fitting and threshold definition

Statistical distribution fitting
A statistical approach is used in this Section.To find the most appropriate distribution for these meaningful defined indices as |neg(LDiñ)|, |LCiñ| and the whole LDiñ distribution, we use the Python stat package, looking for the most appropriate distribution fit.We tested 84 continuous functions implemented in Python (the whole list includes 94: https://docs.scipy.org/doc/scipy-0.18.1/reference/stats.html)plus a powerlaw distribution function package described later.Most of these functions can be found, e.g, in Wilks (2011).The not used distributions due to convergence lack are the following: burr12, exponnorm, gennorm, halfgennorm, invnorm, kappa4, kappa3, levy_stable, skewnorm, trapz, vonmises_line.
These fits and the subsequent intersection points between significant distributions are utilized to define the threshold for intense, very intense and extreme cases.Note that no previous threshold to avoid the distribution bulky part is utilized to fit the functions.All the daily variation usually appears as noise in a histogram.Sometimes the distribution bulk is not regarded for fitting distribution functions; here, since the noise is highly reduced due to the LDiñ processing, the distribution bulk is also suitable for fitting, not only the distribution tail.
The minimization has been done through the Kolmogorov-Smirnoff (KS) distance parameter D computation.The minimization is parametric, through the maximum likelihood estimator (MLE) used to perform the distribution fitting as default; as a double check, we also inspected the minimization by the Kolmogorov-Smirnoff method, producing the same results.
Consequently, the statistics results are robust.Moreover, the fitting is not binning dependent.
Looking for the best fit, we seek the minimum value of the Kolmogorov-Smirnoff distance D for all the distributions for |neg(LDiñ)|, |LCiñ| and the whole LDiñ.Trying to improve the histogram plots, we used the complementary cumulative distribution function (CCDF).Distribution function plotting is made by the survival function computation.The CCDF is appropriate to be fit since it avoids bin dispersion, obtaining an appropriate binning at the distribution tail (as in Riley, 2012).
This bin dispersion appears, for instance, in Fig. 3 at the bottom plots.For now on, the bin size employed from |neg(LDiñ)| is 5 nT and for |LCiñ| is 2 nT min −1 .
For the |neg(LDiñ)| computed CCDF, the best fit is achieved when D is in the order of magnitude of 10 −3 .These distributions are beta-prime (beta of the second kind), the f or Fisher-Snedecor (a type of beta-prime), and the Mielke (beta- For the case of |LCiñ| computed CCDF, the best fit distribution is the power law, obtained with the Python package 'powerlaw' (Alstott et al., 2014).The starting point of the powerlaw In addition to the power law, the Mielke (beta-kappa) is the best fit.Next, the rest of the distributions exhibit a D larger than 0.02, as the Gilbrat (from the lognormal family), Fisk (log-logistic), and the lognormal distribution.The rest of functions appear in Table 3. Unfortunately we did not find any distribution that fit the |LCiñ| with D in the order of magnitude of 10 −3 apart from powerlaw.The best ten fit distributions are displayed in Fig. 5 (middle).
For the whole LDiñ CCDF case, it is more complicated to fit both positive and negative parts.However, the Johnson-Su (Johnson's unbounded, a hyperbolic sine type) is the best distribution, with D of the order of magnitude of 10 −3 .Also the non-central Student's t-distribution (nct) is in the same order.Values of D larger than 0.03 correspond to the Student's tdistribution, the generalized logistic distribution ('genlogistic'), the hyperbolic secant ('hypsecant'), the Laplace, the double Gamma distribution ('dgamma'), double Weibull distribution ('dweibull'), and the logistic.The Cauchy distribution corresponding D is larger than 0.06, and we can see that is departing from the best fitting, as shown in Table 4.The best ten fit distributions are displayed in Fig. 5 (bottom).
For completeness, we mention here the best fit for the whole LCiñ distribution, which is the Johnson Su distribution, with D=0.007; the Student's t has D=0.011; and also the non-central Student's nct gets similar values, D=0.011confirming the LCiñ very small skewness -; the next best fit is the double Weibull, with D=0.026, and Cauchy distribution with D=0.028; therefore general trends are very similar to the whole LDiñ distribution.

Defining thresholds through distribution best fits
As a first approach, we established some arbitrary thresholds in a preliminary way.For |neg(LDiñ)|, these thresholds we used are the percentiles 99%, 99.9% and 99.99%.The corresponding values to these percentiles are -87, -179 and -357 nT, for intense, very intense and extreme, respectively.
In the case of |LCiñ|, setting the thresholds to percentiles 99.9%, 99.99% and 99.999%, the corresponding values are 5, 5 13, 30 nT min −1 .Percentiles are employed, for instance, for the mean residual life of the H derivative for a high latitude observatory, as in e.g., Thomson et al. (2011).However, an advanced way of defining thresholds is making use of the statistical distributions and treatment presented in Section §3.1.For all distributions, we can see that some parts of them have excellent fits and some others whose fitting is not that good, particularly in the distribution tails.Then, taking advantage of this fact, a way of defining the threshold is the intersects of the CCDF with the best fit distribution, or even the intersect of the two best fit distributions.Then, we can define a threshold of a determined degree of geomagnetic disturbance -or its derivative -when the intersect is in the bulk, the tail or the extreme tail.
For the |neg(LDiñ)|, as shown in Fig. 5 (top), the intersects are computed with a relative tolerance of 0.005.Both betaprime (black dashed line) and f (thick red line) intersects the CCDF (in blue) in -100 nT, -205 nT and -475 nT.These values can be set as the thresholds for intense, very intense and extreme disturbances.Lognormal (dashed yellow line) starts departing in -30 nT, and most of the others start departing between -60 and -75 nT, fitting well the bulk but not meeting the heavy tail.This regime difference may mark the threshold setting for moderate disturbances at about -60 nT.
For the case of |LCiñ|, -Fig.5 (middle) -there are two intersects of |LCiñ| shifted CCDF (green line) and the fit powerlaw (in red), at about 6 nT min −1 , which is where powerlaw departs from the shifted CCDF ; and at 18 nT min −1 , where the powerlaw red line meets again the shifted CCDF.We can set the two thresholds in the powerlaw for the case of unsettled and intense values in these two intersects.Next, we can consider the intersect between some the best D-scoring distributions, as Gilbrat (brown line) and betaprime (dashed green line) to determine the threshold for the extreme events.In this case, the intersect is found at 32 nT min −1 , with an absolute tolerance of 10 −10 .The lognormal departs the CCDF at about 6 nT min −1 , clearly separating the bulk from the tail.
The whole LDiñ -Fig.5 (bottom) -is fit for the sake of illustrative purposes on asymmetric distributions, but the threshold in this case is not sought out.

Discussion
An amount of research about distributions applied on different data has been done, such as Dst mainly but also on SYM-H, polar indices, etc. Historically, one of the first applied distributions which may fit the very skewed distribution of geomagnetic disturbances is the lognormal.
Lognormal distribution shapes of the Dst have been commented in Campbell (1996), somehow relating the Dst shape profile (that resembles a lognormal function) to a lognormal distribution.In Aguado et al. (2010), the recovery phase is fit as a hyperbolic function.However, the disturbance profile shapes will not be further discussed here.
Importantly, also quasi-logarithmic scales are used for defining classes, as for the K index (Mayaud, 1980).In Burlaga and Lazarus (2000), distributions of solar wind are fit to lognormal or double-peaked lognormal, since two solar wind distributions are employed, for slow and fast solar wind.Pulkkinen et al. (2012) use shifted lognormal for the geoelectric field amplitudes for past events.Love et al. (2016) employ 1-min horizontal component data from several observatories, fitting to threshold-truncated log-normal with least-squares and maximum likelihood method.Also Love et al. (2015) test both (truncated) lognormal and powerlaw Dst distributions, finding that different statistical tests favour the lognormal.Log-normal and powerlaw distributions are related, since some small changes in generative processes which show log-normal distributions may lead to power-law generative processes (Mitzenmacher, 2003).
Power-law distributions are extensively used in many fields, mainly because they are scale invariant, e.g. in Newman (2005).
A number of power laws appears in X-ray flare energy distribution, solar energetic particle event fluence distribution, auroral activity and radiation belts, due to the consideration of these physical systems as self-critically organized (e.g., Crosby, 2011).
For solar flares, a number of power laws are obtained for flare length, volume, time and diffusivity, revealing the fractal nature and the intrinsic relationship to self-organized criticality (SOC) (Aschwanden et al., 2013).In Yermolaev et al. (2013) the relationship between a variety of interplanetary features and Dst values is studied, fitting log-log quadratic functions, employing power law and square-law fits.Power laws are also utilized for relating the geomagnetic disturbance and its rate of change in scatterplots (Tóth et al., 2014).
Power laws appear as the most utilized distribution in Riley (2012) for CME speeds, X-ray flares and geomagnetic events.
Coincidentally, the Dst distribution index of -3.2 there is on the same order than the power law index for |LCiñ| distribution.
The CME speed distribution also follows a power law with index equal -3.2.Kataoka (2013) fit power laws dH amplitude distribution, finding α around -3.2 (solar cycle 23) to -3.8 (solar cycle 16); this range is in agreement with the results obtained here.Wanliss and Weygand (2007) highlight the lack of characteristic time scale for power laws employed on storm waiting time with SYM-H index and other interplanetary medium parameters.They obtain an index slope of -1.2 using a powerlaw with an exponential cutoff.
Power laws are also the most likely function for power spectra of a myriad of physical magnitudes: for solar wind characteristics, in Burlaga and Lazarus (2000) the power spectra of the solar wind is analysed, finding that the power spectral index for interplanetary medium density is α=-5/3 (coincidental with that of Kolmogorov turbulence index, and also coincidental with the fractral powerlaw for 3D events and volumes (Aschwanden et al., 2013).For proton temperature α also gets values For the |LCiñ|, the powerlaw, Mielke (beta-kappa), Gilbrat (lognormal type), Fisk (log-logistic) best fits correspond to a variety of function families, noting that the fitting is in general worse than |neg(LDiñ)| cases.
For the whole LDiñ, the hyperbolic-sine Johnson-Su has similarity in shape compared to Student's t -and its non-central version nct-, which is also a generalized hyperbolic distribution.
Both f -and beta-prime (inverted beta or beta distribution of the second kind) and power-law distributions are unbounded, in addition to the Johnson-Su (hyperbolic-sine type) and the Mielke (beta-kappa).This also applies to the power-law for |LCiñ|, since it is not bounded either.This means that, statistically, the tails are not limited.Therefore, the best fit functions for |neg(LDiñ)|, LDiñ and |LCiñ| do not have upper limits, theoretically.Then, we may think in the variables plotted in the two-dimensional histogram shown in Fig 4 : the LDiñ and LCiñ may not have any upper limit.This is in agreement with the results in Wintoft et al. (2016).However, the physical parameters that trigger geomagnetic storms may be upper limited by the finiteness on impulsivity of solar features, such as CME speed, solar wind speed, shocks, which cannot take every kind of values, but have an empirical upper limit (finite size effect, e.g.Newman, 2005).
As we mention before, no time between geomagnetic disturbances is considered.As future work, this can be considered for clustering proper geomagnetic storms; then, block-maxima related distributions to those obtained here would be expected.

Conclusions
Seeking the best distribution function to |neg(LDiñ)| and |LCiñ| distributions, we found that the intersect points between the best fit function and the index CCDF can be used to establish thresholds, being more objective than some other ad-hoc thresholds.
The best fits cluster in family distributions: for the |neg(LDiñ)|, beta distribution family is the predominant.In the case of |LCiñ|, there is a variety of families, without a clear predominance.For the whole LDiñ, hyperbolic distributions are the best fit functions.All these functions are unbounded.
(top right), the distribution of LCiñ has a very small skewness.Taking into account the significant values in the form of |neg(LDiñ)| and |LCiñ|, we plotted log-log histograms.In the case of |neg(LDiñ)|, the heavy tail of the distribution is shown clearly in the bottom left of Fig.3.A conspicuous linear trend appears for |LCiñ| histogram in the bottom right panel of Fig.3.We note the bin size variation in logarithmic scale, issue that will be improved in the next Section.We also show a scatterplot of the absolute value of LDiñ and LCiñ to see how their values relate to each other and show the relationship, if any, between extreme values.Unfortunately we were not able to fit a line of this scatterplot due to the large amount of points, similar to those inTóth et al. (2014).Due to the large dataset, we created a two-dimensional histogram instead, plotting the absolute values of LDiñ and LCiñ, in the left panel of Fig.4.The bin size is equivalent to ≈ 6 nT for LDiñ and 0.5 nT min −1 for LCiñ.The histogram evidences that large values of |LCiñ| may happen during non-highly disturbed geomagnetic conditions, according to LDiñ values.An equivalent plot is made in log − log in the right panel of Fig.4, with the same bin amount that the previous one.The truncation in the right plot is at 0.1 nT and 0.1 nT min −1 .Datapoints cluster roughly around 0 for both indices (see red area), which is the noisy data regime.Interestingly, values also clutter in regions (shown in green) where some trends between LDiñ and LCiñ can be noticed.The areas displayed in blue are the most extreme values and they are more scattered.
kappa) distribution.Next, D values are around 0.01 for the Burr distribution (a generalized log-logistic type), exponenciated Weibull ('exponweibull'), generalized Gamma distribution ('gengamma'), powerlaw lognormal ('powerlognorm') and generalized Pareto distribution ('genpareto') and Lomax (a Pareto type-II distribution).Upper D values around 0.02 correspond to the lognormal distribution, and to three other distributions (not shown, and this is marked as [...]).These values are shown in Table 2.The best ten fit distributions are displayed in Fig. 5 (top).
Fig. 5 (middle).The powerlaw use and its low D confirms the first impression on the linear trend of the bottom right panel in Fig. 3.

Figure 3 .
Figure 3. Top left: Normalized distribution of LDiñ values.Bin sizes of 5 nT.Top right: Normalized distribution of LCiñ values.Bin sizes of 0.5 nT min −1 .Bottom left: Log normalized distribution of |neg(LDiñ)| values.Bin sizes of 5 nT.Bottom right: Log normalized distribution of |LCiñ|.Bin sizes of 2 nT min −1 .

Figure 4 .
Figure 4. Left: two-dimensional histogram of the LDiñ vs. its derivative LCiñ.Right: log-log two-dimensional histogram -scatterplot of the the LDiñ vs. its derivative LCiñ.

Table 1 .
This table shows the thresholds of geomagnetic activity for different indices

Table 2 .
D values for best fitting of distributions for the |neg(LDiñ)|.

Table 3 .
D values for best fitting distributions for |LCiñ|.

Table 4 .
D values for best fitting distributions for the whole LDiñ.