On the validation of K index values at Italian geomagnetic observatories

Local K index and the consequent global Kp index are well established three-hour range indices used to characterize the geomagnetic activity. K index is one of the parameters which INTERMAGNET observatories can provide and it’s widely used since several decades, although many other activity indices have been proposed in the meanwhile. The method for determining the K values should be the same for all observatories; so INTERMAGNET consortium recommends 10 a software code, KASM designed for an automatic calculation of K index according to the Adaptive Smoothed method. K values should be independent on the local dynamic response, therefore for their determination each observatory has its own specific scale regulated by the K9 lower limit, which represents the main input parameter for KASM. The determination of an appropriate K9 value for any geomagnetic observatory is then fundamental. In this work we statistically analyze the K values estimated by means of KASM code for the Italian geomagnetic observatories of Duronia (corrected geomagnetic 15 latitude λ~36° N) and Lampedusa (λ~28° N) comparing them with the German observatories of Wingst and Niemegk. Our comparative analysis is finalized to establish the best estimation of the K9 lower limit for these two stations. A comparison of K9 lower limits found for the Italian observatories with results from a previous empirical method was also applied and used to verify the consistency and reliability of our outcomes. 20

As a natural consequence of K index, the planetary geomagnetic activity index Kp was proposed by Bartels (1949). It is derived from the standardized K index (Ks) of 13 magnetic observatories at mid latitude and it is representative of the large spatial-scale of the solar wind-magnetosphere coupling energy. Therefore, K index is the fundamental parameter for Kp estimation; Kp, as any other index, has limitations and drawbacks; however, it's precious since it's an historical parameter and long data series are available; it is widely used, for example, in space weather applications, for identify quietest days (Johnston, 1943) used also in the IGRF modeling, for verifying solar wind driven modulation in the atmospheric parameters during disturbed conditions (Regi et al., 2017).
The main difficulty for K indices evaluation is to assign a proper quasi-logarithmic scale to the geomagnetic fluctuations that satisfy the principle of the assimilation of frequency distributions (AFD): the frequency distributions (or occurrences) of K index values at different sites are, in principle, the same (Bartels et al., 1939). In other words, A values vary from observatory to observatory in such a way that the historical rate of occurrence of certain levels of K is about the same at all observatories (Bartels-Mayaud rules). This implies that, for a given K value, AK increases with increasing latitude, and the fundamental quantity for the K index calculation is represented by the minimum amplitude L9 corresponding to K=9, from which also the other AK values are derived.
The original determination of K indices (Bartels et al., 1939) required hand scaling of analogic magnetograms. For many years, K index was in fact manually scaled by visual determination and removal of the regular daily variation; the remaining largest amplitude of geomagnetic disturbances in the two horizontal components during each 3-hour UT interval, was used to determine the K index values from a conversion table between classes of ranges in nT and K indices.
The question of the derivation of geomagnetic indices from digital data arose at the end of the seventies of the last century.
Different algorithms enabling computer derivation of K indices were then developed and carefully assessed in the frame of an international comparison organized by the IAGA Working Group "Geomagnetic indices" (Coles and Menvielle, 1991;Menvielle, 1991;Menvielle et al., 1995).
This implies the production of computer plots of digital data with scale values similar to those of photographic magnetograms (Menvielle et al., 1995). The International Association of Geomagnetism and Aeronomy (IAGA, http://www.iaga-aiga.org/) promotes tools or methods able to make it easier to keep track of files and analyses done on computers. Different methods were proposed and carefully compared and assessed in occasion of an international meetings organized by the IAGA Working Group "Geomagnetic indices" during the Vienna IUGG general Assembly in 1991 and four methods were acknowledged: FMI (provided by Finnish Meteorological Institute, Finland), LRNS (Hermanus Magnetic Observatory, CISR, South Africa), KASM (Institute of Geophysics, Polish Academy of Science) and USGS (USGS, USA), whose Fortran 77 codes are available at the International Service of Geomagnetic Indices (ISGI, http://isgi.unistra.fr/softwares.php).
We used one of the 4 methods endorsed by IAGA, through the ISGI international service, for calculation of local geomagnetic activity indices K and, in particular, the KASM method that is based on adaptive smoothed method (Nowozyński et al., 1991). For the calculation of the K index, IAGA formatted files are used by KASM code. It requires three daily files, the one under analysis and the files of the previous and following days on which the code estimates the regular daily variation. The code also needs as input parameters the L9 value and the yearly average of the H component relative to the year of interest.
We want to point out that it does not exist a unique L9 at a given geomagnetic latitude since each site maight be affected by different local features such as, for example, crustal anomalies (Chiappini et al., 2000) and/or coast effect (Parkinson, 1962;Regi et al., 2018). Moreover, there is the inevitable MLT dependencies of magnetic disturbances which can be smoothed out trough statistical approach, considering long time observations. For the inclusion of a new geomagnetic observatory into the INTERMAGNET network, an L9 value can be initially assigned according to the ISGI indication but it can be refined by comparing long term geomagnetic field variations at the new observatory and at historical ones for which K indices are estimated by using well defined L9 levels.
We used the geomagnetic data from the two Italian geomagnetic observatories at Duronia (DUR) and Lampedusa (LMP), evaluating the K index with the purpose of estimating the best L9 value for each observatory. Up to now, K index was evaluated only for DUR observatory, using L9=350 nT, both for hand-scaling (since 2012) and for KASM program (since 2017).
In this work we evaluated L9 throughout a correlation analysis performed between K index at DUR with that provided by historical observatories. In order to take into account the magnetic local time dependency of reference K index, European observatories were selected. As possible reference observatories we chose Wingst (WNG) and Niemegk (NGK), since they are among the 13 observatories that contribute to the Kp estimation and their local magnetic time is quite close to that of our Italian observatories.
Our investigations suggest that NGK is the best reference observatory for Italian geomagnetic observatory of DUR, probably due to the closest magnetic local times: indeed the amplitude of magnetic disturbances has dependence on (magnetic) local time which affects the K index values (Chambodut et al. (2013)). By comparing DUR with NGK we estimated a reliable DUR L9 level of 320 nT. Finally, by comparing also LMP with NGK, a reliable LMP L9 level of 310 nT is estimated.

Data and methods of analysis
Geomagnetic field variations at Italian geomagnetic observatories of DUR and LMP are measured by using three-axis fluxgate magnetometers along the northward (H), eastward (D), and vertically downward (Z) directions in the geomagnetic reference frame at 1 s sampling rate. Following the INTERMAGNET directives, geomagnetic time series are also stored as daily archives at 1 min sampling rate, according to the IAGA 2002 format.
In this work we used minute data computed from second data using INTERMAGNET 1s to 1min filter, available in the time interval 1 January 2017 -31 December 2018, a temporal window which falls in the lower part of the sunspot number curve for the cycle 24 (Upton & Hathaway, 2018).
These data are used for estimating K indices by using KASM algorithm which is recommended by INTERMAGNET. In this work, the definitive L9 level at DUR is empirically estimated throughout the following procedure: a) we selected a reference observatory; b) K index time series at DUR are computed by using KASM for different L9 values (KL9) in the range 200-400 nT with a step size of 10 nT; c) each KL9 index time series at DUR is compared with K index time series at reference observatory through correlation analysis; d) the definitive L9 level at DUR is estimated in correspondence of the maximum correlation coefficient.
As possible reference we selected the historical observatories of NGK and WNG since they are among the 13 observatories used for Kp evaluation, they are both in Europe at a MLT close to that of DUR and LMP (see Tab.1 and Fig.1), which is at the moment the principal Italian observatory, member of the INTERMAGNET network and then used in this work as reference observatory for any Italian site. By following our procedure at point c), using independently NGK and WNG, we found that the higher correlation is reached with NGK. The same procedure from a) to d) is applied to the lower latitude Italian observatory of LMP.
We note that K indices at NGK and WNG are both generated by using FMI algorithm. Then, we find useful to verify that FMI and KASM are consistent methods by comparing K values estimated with both methods at NGK.
Finally, we compared L9 values estimated at Italian geomagnetic observatories by means of our method, with those estimated using an historical method proposed by Mayaud (1980).

L9 empirical estimation
We retain that it is important to know how K indices are distributed at consolidate reference observatories of NGK and WNG and how they are in relation to each other. shows also a non-symmetric distribution around zero, with a very different number of cases in correspondence of ΔK=KK=±1: 1103 cases (~19%) for ΔK=KK=-1 and 48 cases (~1%) for ΔK=KK=+1; this feature is also evidenced by the linear regression law: We investigated the frequency distribution ν of K at NGK in correspondence of ΔK=KK=±1 cases, and compared these distributions with the general distribution of K at NGK (the one shown in Fig.2d). Figure 3 shows these distributions, separately for ΔK=KK= -1 and ΔK=KK= +1 conditions during 2017-2018 (top panels), and separately for the two years (middle and bottom panels). In each panel, the K occurrences at NGK, regardless of ΔK=KK, are superimposed (red lines). It can be seen that the general distributions of K and ΔK=KK= -1 are very similar; those of K and ΔK=KK= +1 are also quite similar, although the total number of cases is really lower. Therefore, the occurrences for ΔK=KK≠0 seem not to be led by particular magnetospheric activity conditions. show the smoothed curves computed by using a five-points triangular window, will be used hereafter as actual experimental results for our investigations. It can be seen that the correlation r is higher for DUR-NGK observatories (r~0.915 for L9=320 nT) with respect to DUR-WNG observatories (r~0.908 for L9=290 nT). Regarding LMP, the correlation attains lower values with respect to DUR, with maximum values of r~0.875 for L9=310 nT with NGK, and r~0.870 for L9=300 nT with WNG. The lower correlations between the Italian observatories and WNG could be due to the higher latitudinal separation and possibly to the geomagnetic coast effect at WNG.
As expected, both the L9 limit and r increase with the increasing geomagnetic latitude of referred observatory, here represented by DUR and LMP. In addition, the higher correlations obtained by using NGK are probably due to the lower latitude (i.e. closer to the Italian observatories) and the closer MLT with respect to DUR (table 1). Also at LMP, even if the MLT is closest to that of WNG, the higher correlation is found with NGK: this result suggests that latitudinal effects are dominant with respect to MLT ones. This can be well understood taking into account that the MLT range of all selected observatories is within 11 minutes, well shorter than the 3-hour interval used for K determination.
Anyway, from these results we can assert that for the comparison with Italian observatories NGK is slightly better than WNG, and it will be used hereafter as the reference observatory for DUR and LMP. We can also assume that the best estimation of the L9 value at DUR and LMP is 320 nT and 310 nT, respectively; so the K indices computed by using KASM at DUR with L9=320 nT and at LMP with L9=310 nT represent the best input parameter for the K evaluation for Italian observatories. Figure 5 shows the frequency distribution of the difference between these K-index time series and that computed at NGK.
The occurrences for both DUR and LMP are distributed around zero and in the range [-1:+1]. For a comparison, we show also the difference distribution obtained for DUR using L9=350 nT (dashed blues line); we recall that this is the L9 value we used up to now. It can be seen that, while for L9=320 nT the distribution is almost symmetric around zero, using L9=350 nT it appears more asymmetric, unbalanced towards positive values, confirming that a higher correlation between NGK and DUR is found for L9=320 nT. At LMP the distribution appears slightly asymmetric, and this discrepancy with respect to DUR could be attributed to the larger latitudinal and MLT difference between NGK and LMP.
Since our validation procedure aims to estimate comparable K indices at Italian observatories, we found useful to compute ΔK=KK between DUR and LMP, whose distribution is shown in Fig. 6. It is almost symmetric around zero, closely reflecting the distribution of ΔK=KK computed between K indices at LMP and NGK (from Fig.5); we can also see that |ΔK| ΔK=KK|ΔK| never exceeds 1, confirming the validity of the results obtained at DUR and LMP by using KASM. It should be noted that our comparative investigation is based on K indices at reference observatories of NGK and WNG, which are computed by using FMI algorithm with L9=500 nT, while at DUR and LMP K indices are estimated by using KASM. Therefore, the question arising from our calibration method is: are FMI and KASM algorithms consistent?
In order to answer to this question, we performed a correlation analysis between K indices obtained by using FMI (K FMI) and KASM (KKASM), where the latter K index is obtained for different L9 levels. In this regards we used NGK 1-min data, provided by the INTERMAGNET web site. Since at the moment NGK geomagnetic field measurements are stored as definitive and quasi-definitive data for 2017 and 2018 respectively, we preferred to separate the analysis for the two years. Figure 7a shows the correlation analyses between K indices at NGK (from FMI algorithm) and that computed by KASM by using L9 in the range 350-600 nT, with a step size of 10 nT. It can be seen that, the r maxima (~0.96 and ~0.95) are reached for L9=460 nT in both years. Assuming that KKASM computed for L9=460 nT represents the better K indices in comparison with KFMI, we examined the occurrences of ΔK=KK=KFMI-KKASM (Fig.7, panels b and d). We can see that for ~90% of cases (and for both years) ΔK=KK is equal to zero. Frequency distributions of KKASM (460 nT) and KFMI (500 nT) indices for the years 2017 (panel c) and 2018 (panel e) are also shown. We point out how the distributions are close to each other, suggesting that FMI and KASM are consistent algorithms, even if they are based on different L9 limits applied for the same observatory, in agreement with Coles and Menvielle (1991), Menvielle (1991) and Menvielle et al. (1995).

Comparison with a previous L9 estimation method
As explained in the introduction, the geomagnetic indices are historically assigned throughout visual inspection of magnetograms. The main difficulty for K indices evaluation is to assign a proper value for the L9 limit from which determining the quasi-logarithmic scale to the geomagnetic fluctuations in order to satisfy the AFD principle (Bartels et al., 1939). Mayaud (1968;1980) proposed a method for calculating the geomagnetic activity level L at a given site by comparing the amplitude of geomagnetic fluctuations at the reference observatory (A0) with that, for example, at new one (A) as follows: L=L 0 A / A 0 , where L0 represents the level of geomagnetic activity at the reference observatory, equivalent to L9, and all quantities are dependent on δ= min[ λoval-λ ], i.e. the minimum angular separation between the site, located at geomagnetic latitude λ, and the auroral region, at λoval. According to Mayaud (1968;1980), an approximate value of δ could be given by δ=69°-λ , where the latitude 69° in the corrected geomagnetic coordinate system defines the auroral zone.
We searched a simple relationship which relates L9 (or L) to the geomagnetic latitude of the observatory.
As showed by Mayaud (1980), L9 increases with decreasing δ (L9∝δ -1 ), as expected for a geomagnetic field induced by a current system. Figure 8 shows L9(δ) (blue points) provided in Tab. 5 by Mayaud (1980), considering only northern hemisphere. These points are well represented by a linear law considering an increasing induction effect with increasing parameter x=1/δ. Therefore, by using x=1/δ, previous relationship is linearized and can be formulated as follows L 9 ( x)=αxx+ β . (1) The results of the linear regression analysis performed on the experimental points are also reported in Fig.8. Equation 1 is therefore useful for estimating a reasonable L9 limit at a different site. In order to evaluate L9 at DUR, LMP and, for a comparison, at NGK observatories, the corresponding δ parameter is required. However, it is not clear how δ was estimated by Mayaud (1980), since it requires, for example, an auroral oval model for estimating the λoval, and an IGRF model for evaluating the geomagnetic latitude λ of a given site (this aspect will be further discussed at the end of this section).
In this regard, we empirically estimated δ(λ) by a linear fit of the experimental data reported by Mayaud (1980). Figure 9 shows experimental points (red stars) and the corresponding linear law (red line) which allows us to extrapolate an estimation of the theoretical δth(λ) for the observatories of DUR, LMP and, for comparison at NGK too (blue circles), where λ represents the corrected geomagnetic latitude used by Mayaud (1980). Finally, by inserting δth into the Eq. (1) we estimated the L9(δth) level at the observatories of DUR, LMP and NGK.
All these results are reported in Tab.2, which also shows for a comparison, L9 obtained by computing δa=69°-λ, and L9exp experimentally derived by our calibration procedure (L9exp), together with the 95% confidence intervals for the fitted L9 values.
It can be seen that all L9exp are consistent to each other within their respective confidence interval, at a given observatory.
The small difference between L9exp and L9th could be due to a different method for calculating the geomagnetic coordinates used in Mayaud and in this work (we use AACGM). In order to verify this hypothesis, we performed a correction on the key parameter δ(λ) as follows: we computed the AACGM latitudes Λ of geomagnetic observatories from Tab. 5 of Mayaud and corrected λ through linear relationship λC= l Λ + m; we performed a linear fit of δ(λ), λC, which provides the relationship for the adjusted δA(λC); finally we performed a linear fit of L9(δ), δA which provides the adjusted L9A(Λ) estimated at our geomagnetic observatories as a function of AACGM latitude.
With respect to the L9(δ) it can be seen that the adjusted L9A(Λ) (shown in Tab.2) are closer to the experimental L9exp, indicating that the correction on geomagnetic coordinate makes a significant contribution on the L9 estimation. LMP is the only one that shows a discrepancy between L9exp and L9 here estimated with different methods. A possible reason of this discrepancy lies in the low latitude of LMP observatory where the ring current and/or electrojet currents dynamics could affect L9 estimations.

Discussion and Conclusions
Four different automated methods for calculation of local geomagnetic activity indices K were endorsed by IAGA, through the ISGI international service, and distributed from ISGI (http://isgi.unistra.fr/softwares.php) web site. For the Italian geomagnetic observatories of Duronia (DUR) and Lampedusa (LMP) we used one of them, i.e. the KASM algorithm that is based on adaptive smoothed method (Nowozyński et al., 1991) which is the only one provided also by INTERMAGNET (https://www.intermagnet.org/publication-software/software-eng.php) An input parameter required by KASM code, as well as FMI code, is the L9 value, the so-called "K=9 lower limit", which allows to determine, for each magnetic observatory, the conversion table between classes of geomagnetic field variation ranges and K index values.
We found L9 values for DUR and LMP through a correlation analysis using as reference the corresponding data from the two European observatories of Wingst (WNG) and Niemegk (NGK), both located in Germany. The choice of these two observatories was prompted by the fact that they are among the 13 observatories which provide their K indices for the determination of the planetary Kp index and moreover, their magnetic local time is very close to that of the Italian observatories.
We note that NGK is the best reference observatory for Italian geomagnetic observatory of DUR, possibly due to the closest magnetic local time; indeed the amplitude of magnetic disturbances has dependence on (magnetic) local time which inevitably reflects on different K index values (Chambodut et al. (2013)). regardless if manually or automatically computed. Our analysis also highlighted the possibility of establishing a linear relationship between a pair of analyzed observatory datasets that can be useful for predicting or deriving the index of one when the other is known.
Another interesting result that we found is related to the consistency of the KASM code and the FMI code, the latter in use at the two German observatories for the K index computation and subsequent release. Although FMI code is based on a different procedure, we verified that the results obtained are consistent with those obtained by KASM code and stable in the two-year time interval, although with a slightly different value of the input L9 parameter. This confirms that the choice of a certain algorithm in place of another does not invalidate the results.
Before the introduction of automatic procedures, based on the definition introduced by Bartels et al. (1939) for the K index concept, in the '80s of the last century Mayaud (1980) used an empirical relation to calculate the level of the local magnetic activity L (equivalent to the L9 values) for a generic point of observation with respect to a referenced observatory. Through a linearization process, we used this relation, which includes some approximations and the necessity of determining the minimum angular separation between the observational point and the auroral region, i.e. a method for determining the geomagnetic latitude, obtaining an independent estimate of the L9 values for our observatories which is consistent, within the 95% interval of confidence, with that obtained by our previous analysis. Moreover, Mayaud (1980) notes that the limitation of the method he proposes is that it is conceived for sub-auroral and mid latitudes; indeed, he suggests that for lower latitudes a constant L9=300 nT can be chosen. This very approximate value is not very far from the values we estimate (320 nT for DUR and 310 nT for LMP), but would certainly be not accurate as them in the comparison with the values from other reference observatories: indeed our results clearly show that a very precise L9 limit is necessary for obtaining K values well consistent at different sites. As a final remark, from the overall view of this work, we are also convinced that the habit to round the value of L9 in multiples of 50 nT, firstly suggested by Bartels et al. (1939) cannot be changed for observatories which have been providing the K index in the past since homogeneity of the series is of primary importance. However, before a new observatory starts providing K index values, it is worth evaluating K indices with a provisional L9 value, assigned by ISGI, and then refine with a procedure like the one we have shown, based on the comparison with reference, well set, observatories.

Author contribution
DDM, SL and MR planned the study. MR performed the data analysis and wrote the codes. PB studied the KASM code functionality and, together with MR, tested and validated its results. DDM, LC and SL improved the quality of the manuscript. All authors read and approved the final manuscript.     20 Table 1. Geomagnetic observatories used in this study; geographic coordinates; Altitude Adjusted Corrected Geomagnetic Coordinates (AACGM), estimated by using Shepherd (2014) algorithm at 100 km above the observatory and magnetic local time at 0 UT.  Table 2. L9 estimated by different procedures. δa=69°-λ (approximate angular distance from auroral region); L9(δa) obtained from δa using linear fit in Fig.8; δth (δ estimated from AACGM lat and the linear fit in Fig. 9); L9(δth) obtained from δth using linear fit in Fig.8; L9A (Λ) obtained from correction procedure on the key parameter δ(λ) explained at the end of section 3.2. The 95% confidence intervals are also indicated.