The Geomagnetic Kp Index and Derived Indices of Geomagnetic Activity

The geomagnetic Kp index is one of the most extensively used indices of geomagnetic activity, both for scientific and operational purposes. This article reviews the properties of the Kp index and provides a reference for users of the Kp index and associated data products as derived and distributed by the GFZ German Research Centre for Geosciences. The near real‐time production of the nowcast Kp index is of particular interest for space weather services and here we describe and evaluate its current setup.

In space weather application, Kp is used for the National Oceanic and Atmospheric Administration (NOAA) Geomagnetic Storms or G-scale, ranging from minor (G1) for Kp = 5 to extreme (G5) for Kp = 9 (www. swpc.noaa.gov). It is also one of the parameters distributed by ESA's Space Situational Awareness Program (http://swe.ssa.esa.int/). Another example for space weather application is the radiation belt forecast system described by Horne et al. (2013), which is based on Kp predictions. It covers the outer radiation belt and the slot region between the inner and outer radiation belt and gives radiation risk indices for internal charging of satellites in geostationary and medium Earth orbit. The Kp index is also used to describe and predict scintillations in global navigational system signals at high latitudes. These scintillations are due to strong plasma gradients of ionospheric plasma irregularities (e.g., Xiong et al., 2016Xiong et al., , 2018. Secan et al. (1997) used Kp to describe and forecast the intensity, latitude and local time of irregularities with the Wideband Ionospheric Scintillation Model and Kp has high skill and impact to machine learning-based predictions of high-latitude scintillations (e.g., McGranaghan et al., 2018). Geomagnetic activity affects thermospheric density (e.g., Bruinsma, 2015) and Kp plays an outstanding role in low earth orbit space safety and space traffic management (Berger et al., 2020;He et al., 2020). Finally, Boteler (2001) showed a relation between Kp and modeled peak electric field magnitudes due to geomagnetically induced currents (GIC).
In more than 70 years since its inception, Kp has proven to be an important and reliable index. Still, some limitations exist compared to other, younger geomagnetic indices. For example, the time resolution of Kp is lower than for the PC (Polar Cap) and AE (Auroral Electrojet) indices (one minute resolution) or the Dst index with one-hour resolution (e.g., Rangarajan, 1987). Also, Kp is capped to a certain value (i.e., 9) for events with extreme geomagnetic disturbance. In this regard, new Kp-like but open-ended indices with 30 and 60 min time resolution have been developed as part of the European Union (EU) project Space Weather Atmosphere Models and Indices (SWAMI; Jackson et al., 2020). These indices, which are already available on https://www.gfz-potsdam.de/en/hpo-index/, will be described in a separate paper. Also, Kp does not reflect the universal time (UT) dependent variation of geomagnetic disturbance , while the Kp-like am index contains this information (Mayaud, 1980;Svalgaard, 1976).
A nowcast Kp is crucial for space weather applications as described above and for space weather monitoring as it warns operational users in near real-time of potential space weather impacts. GFZ German Research Centre for Geosciences (GFZ) has a nowcast system, which will be later described and discussed in comparison to the definitive Kp index and the nowcast system by NOAA.
The Kp index is available back to 1932 and thus valuable for space climate studies. With the present publication serving as a reference, the Kp index and derived products are now distributed with a digital object identifier (DOI; Matzka et al., 2021, data publication) and under the CC BY 4.0 license.
In what follows, we give a first hand account of the current derivation and distribution of the geomagnetic Kp index and its derived products (sections 2.1 and 2.3) and its historical development (section 2.4). We demonstrate some of the fundamental properties of Kp based on the data from 1932 to 2019 in section 3.1 and investigate the influence of the Kp-stations' change in geomagnetic latitude due to secular variation (section 3.2). In section 3.3 we describe the production and properties of the nowcast Kp index. Section 4 describes the distribution and license of the index and section 5 gives a summary and conclusion. We mostly use the terms and notation of Bartels (1957aBartels ( , 1957b. Like him, we use italics for the letters denoting indices, however, we skip the hyphen in the index names (we use, e.g., "K index" and "Kp index").

Derivation of K
Note that K indices are a standard product of many geomagnetic observatories from all latitudes, while only 13 subauroral observatories contribute with their K indices to Kp. More details on geomagnetic observatories are given in, e.g., Matzka et al., (2010) or Chulliat et al. (2017). K is an observatory specific geomagnetic index defined on the eight three-hourly intervals of the UT day 00-03, 03-06, …, 21-24 UT . K is determined from the two horizontal geomagnetic components. Until 1963, also the vertical component Z has been used to derive K (see section 2.4 for the history of K and Kp and changes in their production). Observatories have different configurations to orient the two perpendicular sensors that measure the horizontal geomagnetic variation. With digital data, the geomagnetic field variations along any field component can be easily calculated from any of these configurations, but in the times when photographic recordings were used, typically until the 1980s, the K indices had to be determined directly from photographic paper and thus from the field components that were actually recorded. In the XYZ-configuration, the horizontal geomagnetic field variations are recorded in the geomagnetic component X along geographic north and Y along geographic east. The other typical configuration is the HDZ-configuration. Here, the horizontal geomagnetic variations are recorded by sensors that are mounted approximately along local magnetic north (as an approximation for the variations of the horizontal field strength, H), and local magnetic east (as an approximation for the variations, in nT, of the geomagnetic declination).
K is derived by first estimating the quiet curve (black curve in Figures 1a and 1c) for each of the two horizontal geomagnetic components (blue in Figures 1a and 1c) at a geomagnetic observatory. In the absence of an SFE (regarding SFEs, see section 2.4.2), the quiet curve is a good representation of the non-K-variation and subtracting it from the geomagnetic variations yields the K-variation (blue curve in Figures 1b and 1d). Various algorithms exist to estimate the quiet curve and individual observatories made different choices regarding algorithms (see section 2.4.3 for more details on the determination of the quiet curve). The K-variation's range r (in nT, that is the same that was labeled "a" by Bartels), i.e., the absolute value of the difference between the K-variation's maximum and minimum, is determined separately for each three-hourly interval and for each horizontal component. In Figures 1b and 1d the red box indicates the maximum and minimum for each interval and the value of range r is given by the red number at each red box. Additionally, a double-headed vertical arrow indicates the range r in the interval 21-24 UT.
In each three-hour interval, the K index is determined by the horizontal component with the larger range r by converting r to K (panel e in Figure 1) according to an observatory-specific scale of lower range-limits, also called the range-limits, and denoted here as K-limits, or individually as the K0-(for K = 0) to K9-limit (for K = 9). Table 1 lists the K-limits chosen for Niemegk (Bartels, 1938), which is regarded as the standard observatory for K indices . For example, the range r with 40 nT ≤ r < 70 nT will be assigned K = 4 and the K9-limit is 500 nT. Additionally, Table 1   (a) X component geomagnetic variation (blue) and non-K-variation (quiet curve, black). (b) X component K-variation (blue), a red box indicates its maximum and minimum in each three-hour interval. The red numbers on top of each box indicate the corresponding range r in nT, which is also indicated by a vertical arrow for 21:00-24:00 UT. (c) and (d) same as (a) and (b), but for the Y-component. (e) K indices for Niemegk for the respective three-hour intervals, as derived from the larger of the corresponding ranges r in panel (b) and (d) according to Table 1.  were pleased by the K-limit scale chosen for Niemegk because the intermediate values K = 2, 3, and 4 occur in the majority of cases (>50%) and both lowest (K = 0 and 1) and higher degrees of disturbance (K > 4) are distinguished well. This frequency distribution is still well maintained for K at Niemegk and for Kp for 1932-2019, see Table 1.
By definition, the K-limit scales for all observatories are proportional to the Niemegk scale (Table 1). The scales are referred to as quasi-logarithmic (Bartels, 1957b). They start with the K0-limit of 0 nT, then K-limits are doubling from K = 1 up to K = 4 in the logarithmic part of the scale, and increase less steeply for K > 4. The K1-limit is a hundredth of the K9-limit.  determined the individual K9-limits of other observatories such that their frequency distribution of K values is closest to that of Niemegk. Practically, this was achieved by choosing the K5-limit such that the frequency of K ≥ 5 for any observatory is closest to that of Niemegk, which is slightly more than 5% (Table 1). According to Bartels (1957aBartels ( , 1957b, only scales with certain K9-limits, namely 300 (for equatorial observatories), 350, 500, 600, 750, 1,000, 1,200, 1,500, and 2,000 nT, should be used, but not all observatories adhered to this rule and some adopted scales based on other K9-limits. In general, the observatory specific scales are strongly dependent on corrected geomagnetic latitude, or more specifically, on the latitudinal distance to the auroral zone (Mayaud, 1980) and from the 1960s Mayaud used the relationship between K9-limit and auroral distance to assign K9-limits to geomagnetic observatories.
The derivation scheme for K comes with some degree of ambiguity and approximation, e.g., by using geomagnetic recordings in the XYZ-configuration or HDZ-configuration, by using photographic recordings or digital data of various time resolutions, by the estimation of the quiet curve, or by the choice for pre-defined discrete levels of K9-limits. This requires extra precaution at a geomagnetic observatory to maintain the reproducibility and homogeneity of its K index time series. Any change in procedure at a specific observatory needs to be checked for consistency of results, typically by employing the old and the new procedure in parallel for a certain period of time. Certain aspects of changes are discussed in section 2.4.
The three-hourly ranges and consequently the K indices are local time dependent (with a maximum in geomagnetic disturbance around local midnight) and thus the mean frequency distribution of K indices is different for the eight intervals of the day (Bartels, 1949;Bartels et al., 1940). The frequency distribution of K values is also station dependent.

Definitive Kp, ap and Ap
Kp is currently derived from 13 contributing geomagnetic observatories, the Kp-stations shown in the upper panels of Figure 2. Table 2 gives the current Kp-stations together with their operating institutes, the year when they became Kp-station, their geographic longitudes and quasi-dipole (QD) latitudes, their K9-limits and respective weight in Kp. QD coordinates were calculated using IGRF-13 (Alken et al., 2021) and the Fortran code provided by . The results depend on the maximum degree and order L of the spherical harmonic expansion, and the maximum degree N of the vertical polynomial expansion. Here, we have used L = 6 and N = 8. The QD coordinates Laundal & Richmond, 2017) for the Kp-stations are almost identical to corrected geomagnetic (CGM) coordinates (Gustafsson et al., 1992) used in many earlier studies.  . b from . c Frequency single-digit Kp, see section 2.2.

Table 1 K-Limits and Frequency Distribution of K for Niemegk and of Kp
Like others (e.g., Kauristie et al., 2017;Takahashi et al., 2001;Thomsen, 2004), we regard all Kp-stations as subauroral. We note, however, that the subauroral latitudinal range is not clearly defined and depends on the phenomena under investigation and various ranges have been given from, e.g., 40°-55° CGM latitude (Mayaud, 1980) to, e.g., 60°-65° dipole latitude (Oguti, 1993). The choice of subauroral latitudes for the Kp-stations is not only crucial for the scientific meaning of the index, e.g., as a proxy for magnetospheric convection (Thomsen, 2004), it also comes with technical advantages for the index production. The separation of K-and non-K-variation is clearer at subauroral than at lower latitudes because, first, the auroral K-variations are more distinct, second, the quiet curve is easier to predict (for the northern component it MATZKA ET AL. 10.1029/2020SW002641 6 of 21  has less day-to-day variability since it lies further away from the mid-latitude Sq foci), and third, during the main phase of a geomagnetic storm, the ionospheric auroral K-variations dominate the ring current signal, which is a non-K-variation. Using measurements at subauroral latitudes also has advantages over the data from auroral latitudes. At auroral stations, the K-variation is mostly of local ionospheric polar electrojet origin (and predominantly in the northern horizontal component, which then determines the K index). At subauroral stations, K-variation from large scale ionosphere-magnetosphere current systems (e.g., McPherron & Chu, 2017) play an important role (and show up, with comparable range, in both the northern and eastern horizontal component). Therefore, K indices from subauroral stations are sensitive to a much larger current system and broader longitudinal sector of the auroral oval than K indices from auroral stations and hence depend less on regional scale geomagnetic activity, see p. 29 in Mayaud (1980). Consequently, it is expected that the global geomagnetic activity level can be gauged by a smaller number of subauroral stations than would be required from auroral stations. At the same time, the subauroral latitudes are more suitable to magnetically monitor the equatorward auroral oval expansion during strong storms as the auroral oval moves away from the auroral stations, which might see a concurrent decrease in K-variation.
The definitive Kp index is derived from the K indices of the Kp-stations, which are sent to the Adolf-Schmidt-Observatory for Geomagnetism Niemegk on a semimonthly basis. There, the K indices are converted to standardized Ks indices by using conversion tables. Following Bartels (1949), in these tables, integer values from 0 to 27 are given for each possible value of K for each UT three-hour interval and Lloyd season (see Table S1, including a short description of their use). The three Lloyd seasons are equinoxes (March, April, September, October (MASO)), June solstice (May, June, July, August (MJJA)) and December solstice (November, December, January, February (JFND)). The integer values 0 to 27 represent 3*Ks. Dividing 3*Ks by 3 gives Ks, an index expressed on a scale of thirds ranging from 0 to 9. Bartels (1949) determined the conversion tables such that the frequency distribution of the Ks values for each Kp-station and in each threehour interval of the UT day best resembles the so-called standardizing distribution. This standardizing distribution corresponds to the frequency distribution, averaged over all Kp-stations, of the K indices in the two UT intervals closest to local midnight and in the respective Lloyd season. Deviating from this principle, K = 0 is always assigned 3*Ks = 0 and K = 9 is always assigned 3*Ks = 27. Conversion tables exist for each Lloyd season (JFND, MASO, MJJA) to account for the seasonal change of the local time dependence of K. At the season boundaries, the seasonal change is smoothed with the help of transition tables (   Bartels, 1949). Additionally, K = 0 around midnight was found to be less frequent in summer (MJJA in the northern hemisphere, JFND in the southern hemisphere) than in winter. Combined with the larger number of Kp-stations in the northern hemisphere, this would lead to an artificial seasonal variation. This artificial seasonal variation is compensated in Kp (Bartels, 1949) through the conversion tables for MJJA (note that for K = 1 and 2 the values for 3*Ks are smaller for MJJA than for JFND). Also, on purpose, no 3*Ks value appears twice in the conversion tables within any eighth of the UT-day for any observatory, allowing K to be reconstructed unambiguously from Ks. The average 3*Ks of Eyrewell and Canberra as well as of Uppsala and Brorfelde count as one Kp-station each (weight in Table 2). The average of 3*Ks for a specific three-hour interval is rounded to the integer 3*Kp, i.e., three times the Kp index. Figure 2 (lower panels) shows, as an example, the Kp index for the geomagnetic storm on September 7 and 8, 2017, and preceding days.
The definitive Kp index as described in the previous paragraph is endorsed by the International Association of Geomagnetism and Aeronomy (IAGA). It is available back to 1932. Kp is a global index of geomagnetic activity for each three-hour interval of the UT day. It is traditionally represented by either 3*Kp = 0, 1, …,27, or in symbol notation by Kp = 0o, 0+, 1−, 1o, 1+, …, 9−, 9, where o, + and − represent the integer (o), plus one third (+) and minus one third (−), respectively. Sometimes, Kp is represented as Kp = 0, 0.3, 0.7, 1.0, 1.3, …, 8.7, 9.0 (or without decimal points: 0, 03, 07, 10, 13, …, 87, 90), and one must not misinterpret these as direct numerical values (e.g., 0.3 is not identical to one third). A traditional graphical way to represent Kp is the so-called note-script or musical diagram (Bartels, 1949), which is organized according to Bartels solar rotation and is available back to 1932 (see section 4. and an example in Figure S1). Kp may sometimes also be rounded to a single-digit integer number and presented with reduced resolution as 0, 1, …9. Examples for the use of the single-digit Kp are the comparison of the frequency distribution of Kp with that of K at Niemegk in Table 1.
Kp with its underlying quasi-logarithmic scale does not lend itself for the calculation of arithmetic means (e.g., daily means, yearly means). To this end, Kp is converted to the linear ap, the three-hourly equivalent planetary amplitude according to Table 3. Sometimes, and to a certain degree confusingly, ap is assigned a unit of 2 nT, sometimes it is regarded as an index without units (Siebert, 1971). By multiplying the numerical value of ap by 2 nT, the average range in nT of the geomagnetic disturbance at a Kp-standard station (located, like Niemegk, at about 50° geomagnetic latitude, with K9-limit = 500 nT) is obtained. For example, MATZKA ET AL.  Bartels et al. (1957b). b Siebert (1971) recommends to regard ap as dimension-less index and to ignore the unit of 2 nT. c Frequency distribution 1932-2019. we have ap = 80 for Kp = 6 (see Table 3) and 80*2 nT = 160 nT is the mean disturbance range for K = 6 at Niemegk (120 nT ≤ r < 200 nT, see Table 1). The average (arithmetic mean) of the eight ap indices of a UT day rounded to an integer is the daily equivalent planetary amplitude Ap (again in units of 2 nT, or dimensionless). Monthly and yearly means of ap have also been labeled "monthly Ap" and "yearly Ap".

Derived Products Cp, C9, and International Quiet and Disturbed Days
The planetary international magnetic character figure Ci is an index predating Kp and not in use anymore. Based on the daily sum of ap, which has a higher numerical resolution than the rounded Ap values, for each UT day, a planetary geomagnetic index Cp = 0.0, 0.1, … to 2.5 is derived that replaces the earlier Ci index (for details, see Mayaud, 1980;Rangarajan, 1987). Cp is still a standard product and we publish it along Kp in some data products. The Cp index can be contracted to the single digit integer index C9. Tables to convert the daily sum of ap to Cp and C9 are given by Bartels (1957b) and in Tables S2 and S3. The so-called R9C9-plots ( Figure 4 in Siebert & Meyer, 1996) showed both C9 and the three-day mean sunspot number and are now replaced by SN9C9-plots to reflect the change in the sunspot number series in 2015 (Clette and Lefèvre, 2016), see section 4). The SN9C9-plots allow direct graphical comparison of sunspot number and C9, ordered according to solar rotations (also available back to 1932, see section 4. and an example in Figure S2).
The International Quiet and Disturbed Days are determined from Kp on a monthly basis. For each month, days are ordered according to the daily Kp sum, the daily sum of the square of Kp, and the largest Kp of the day (Bartels, 1957b). The days in the ordered sequences are then numbered and for each day the mean of the three numbers is a measure of its relative geomagnetic activity within this month. The month's 10 most quiet days and the 5 most disturbed days according to this mean are published by GFZ as the International Quiet and Disturbed Days (see section 4). They are available back to 1932. The list of International Quiet and Disturbed Days is often applied to identify geomagnetically quiet or disturbed periods for geomagnetic and space physics studies.

Historical Development of Kp
The previous sections describe the current definitive Kp index production. Here, we describe historical changes to the Kp index and to the rules for its production.

Associations and Institutions Involved
The K index was adopted in 1939 by the International Association for Terrestrial Magnetism and Electricity (IATME), the predecessor of IAGA.

Changes of Rules Affecting K and Kp
Until 1963, the K indices were determined from the geomagnetic disturbance in all three components. However, geomagnetic induction (e.g., Kuvshinov, 2008) enhances the variation of the horizontal components and decreases the variation of the vertical component Z. And since only the component with the largest range determines the K index, the Z component is not expected to play a role for K except for stations in close proximity to the auroral electrojet and for stations with anomalous subsurface induction effects (Siebert, 1971). Induction also smoothens the variation in the Z-component, making it more difficult to separate K-variation from non-K-variation in Z (Mayaud, 1980). IAGA resolution 4 from 1963 determined "…that from the 1st of January 1964, the Z-component will not be used for the measure of the three-hourly K-index, except by the standard Kp-observatories, …". The advent of computer-derived Kp indices and its implication for the Kp-stations in this context will be discussed in section 2.4.3.
SFEs can be difficult to identify in magnetograms even for experienced observers and hence, from 1968 onwards (according to IAGA resolution 12 in 1967), the K and Kp indices should include SFE's as K-variation instead of treating them as a non-K-variation. This simplified the production of K and Kp (also for the later introduced computer algorithms, see section 2.4.3), and has little practical influence on the Kp series due to the very rare occurrence of SFEs. For time intervals affected by an SFE, an additional K′ and Kp' index should be determined that treats SFE's as non-K-variation. K′ and Kp' should be published supplementary to Kp. Note that the same instruction was already published 10 years earlier by Bartels (1957a). Note also that there currently is a lack of procedures for collecting K′ as well as for distributing and archiving Kp'. This should be addressed in the future.

Introduction of Algorithm-Derived K Indices
The most complex task in deriving K indices is estimating the non-K-variation (or quiet curve) in order to determine its counter-part, the K-variation (or geomagnetic disturbance). Since Sq has a significant day-today variability, estimating the quiet curve for hand-scaled K indices required experienced personnel familiar with the typical geomagnetic variation at a station. More than 95% of the K indices obtained independently by two persons would typically be identical and any differences in the derived K indices would be randomly distributed and never exceed one unit (Menvielle, 1981;Rangarajan, 1987). Experienced persons that are hand-scaling K indices for observatories with which they are not personally familiar would achieve 92% identical results when excluding ranges r close to a K-limit, and 80%-85% including the cases with ranges r close to a K-limit. It was argued that the difference between algorithm-derived K indices and hand-scaled K indices should not be larger than the typical difference between independently hand-scaled K indices. Based on one year of data from 10 subauroral and 1 auroral observatory, IAGA finally endorsed four algorithms that achieved at least 69.9% accordance to hand-scaled K indices, with differences in derived K indices greater than one amounting to less than 2% (Menvielle et al., 1995). The four endorsed algorithms are the Adaptive Smoothing Method (ASM; Nowozynski et al., 1991), Finnish Meteorological Institute method (FMI; Sucksdorff et al., 1991), linear-phase robust non-linear smoothing (LRNS; Hattingh et al., 1989) and United States Geological Survey method (USGS; Wilson, 1987). The algorithms in use today at the various Kp-stations are FMI (for observatories NGK, UPS, WNG), semiautomatic LRNS (for CNB), USGS (for FRD, SIT), two modified ASM algorithms of which one is BGS (for HAD, ESK, LER) described in Clark (1992) and the other one is used for MEA and OTT, and additionally the NZ algorithm (McNoe, 1989) for EYR. Geomagnetic observatories traditionally are very concerned about the homogeneity of their data series and since all these organizations are aware of the importance of their K indices, we can assume that the agreement between algorithm-derived and hand-scaled K indices was carefully tested at each Kp-station. This test is documented in Wilson (1987) for the USGS method and Clark (1992) for the BGS method. The algorithms use one-minute observatory data, typically from several days, which is smoothed to yield the quiet (or smooth) curve for a particular UT day. However, all these algorithms have in common that K is determined from the horizontal components alone, and not from the Z component, even for the Kp-station closest to the auroral zone. The provision made in IAGA resolution 4 from 1963 that Kp-station should continue to produce their K indices based on all three components is not any more adhered to with algorithm-derived K indices. The effect on K and consequently Kp should nevertheless be minor (see Section 2.4.2).

Changes to the Network
Some Kp-stations had to be replaced for operational reasons.  (Siebert & Meyer, 1996). In 2004, Uppsala replaced Lovö as Kp-station in Sweden. The current Kp-stations are listed in our Table 2 and the current and previous stations are shown in the upper panels of Figure 2.
While Agincourt had a K9-limit of 600 nT, Ottawa has a K9-limit of 750 nT due to its higher QD latitude. Note that Canberra, which replaced Toolangi (K9-limit 500) in 1981, has a K9-limit of 450 nT, a value not originally foreseen as K9-limit by Bartels (1957aBartels ( , 1957b. All other replacing (new) Kp-stations have the same K9-limit as their predecessor station. When Kp-stations were replaced, no changes in the calculation of Ks from K were introduced to compensate for the change in location (Rangarajan, 1987), except for Niemegk replacing Witteveen (see Table S1).
The magnetic field of single-wire, high-voltage DC power lines between Germany and Sweden affect Brorfelde (BFE) observatory since December 1994 ). The effect on K was investigated from December 2008 to February 2009, a period of very low geomagnetic activity with 30% of the K values of BFE being 0 (mean Ap < 5, corresponding to Kp < 2-). Errors in K derivation due to the power line magnetic field predominantly had the effect of shifting K = 0 to K = 1. This happened for 13% of the K = 0 values or 4% of all K values (Fox- Maule et al., 2009). For higher magnetic activity, the rate of K = 0 and hence the rate of affected K values will decrease. Moreover, since Brorfelde contributes with weight 0.5 to Kp, the overall effect on Kp is expected to be small. Table 2 and the effect of changes to the network will be further addressed in section 3.2.

Fundamental Properties of Kp
This section describes basic characteristics of Kp and its relation to other solar-terrestrial parameters, which help users understand the temporal behavior of Kp. The long-term variation of Kp is illustrated in Figure 3a as a time series of annual mean Kp values from 1932 to 2019. The annual means were calculated from the linearly scaled ap index and then converted back to the equivalent Kp values. Also shown in Figure 3a is the total sunspot number for the corresponding period. The sunspot number exhibits the well-known solar cycle with a period of ∼11 years. The annual mean Kp also shows an 11-year period, but tends to be largest during the declining phase of the solar cycle, i.e., when overall geomagnetic activity is elevated due to long-lasting high speed solar wind streams (Lockwood et al., 1999). It can be noted in Figure 3a that geomagnetic activity was particularly low during the solar minimum year of 2009, which resulted from unusually weak solar wind forcing (Russell et al., 2010).
Kp is strongly connected with some solar wind parameters. Figure 3b shows the dependence of Kp on solar wind speed and interplanetary magnetic field (IMF) clock angle. The IMF clock angle is defined as θ = arctan (By/Bz), where By and Bz are the eastward and northward components of the IMF, respectively. For θ = 0˚, IMF is purely northward in the Y-Z plane in geocentric solar magnetic (GSM) coordinates. Similarly, IMF is purely eastward for θ = 90˚, westward for θ = −90˚, and southward for θ = ±180˚. We used the OMNI 1-min solar wind data estimated for the Earth's bow shock nose (King & Papitashvili, 2005) for the years 1981-2019. The data were time-shifted by 20 min to take into account the time for solar wind signatures to propagate from the bow shock nose to the ionosphere (Manoj et al., 2008) before being averaged over three-hour intervals and compared with the corresponding Kp. As seen in Figure 3b, high Kp occurs preferably under southward IMF conditions (θ = ±180˚). This is because the southward IMF can most effectively merge with the Earth's northward magnetic field, which leads to solar wind energy injection into the magnetosphere. As shown in Figure 3b, Kp also depends on the solar wind speed, as the merging rate increases with increasing solar wind speed. A comprehensive discussion was presented by Newell et al. (2007) and Kauristie et al. (2017) as to how Kp and other geomagnetic indices are related to various solar wind parameters.
Kp variations on intra-seasonal time scales are dominated by the solar rotation effect, which produces a spectral peak around 27 days, as seen on the Lomb-Scargle periodogram in Figure 3c. Depending on the zonal distribution of active regions of the Sun, spectral peaks also appear at ∼13.5 and ∼9 days. These periodic variations have a significant influence on the weather of the thermosphere and ionosphere e.g., Lei, Thayer, Forbes, Sutton, & Nerem, (2008); Lei, Thayer, Forbes, Wu, et al. (2008).
Local geomagnetic activity and hence K indices are known to depend on local time (LT) and season. Averages of K indices from a globally well-distributed network of stations should exclude local time influences  but might depend on universal time UT due to processes like the daily rotation of the Earth's geomagnetic dipole axis with respect to the Sun-Earth line. The standardization process described in section 2.2 should remove, to a large extent, the UT dependence of Kp . In fact, the dependence of Kp on UT is only one third of a unit (during equinoxes) or smaller (Figure 3d). Kp, however, shows a marked seasonal variation with equinoctial maxima (Figure 3d) with an amplitude of about one unit. This is probably due to geomagnetic storms that are more frequent in equinoxes than in solstices (e.g., Cliver et al., 2000;Russell & McPherron, 1973).
Thus, Kp contains information both on the solar wind as well as on the coupling between the solar wind and the magnetosphere-ionosphere system (Cliver et al., 2000;Russell & McPherron, 1973;Siebert & Meyer, 1996).

Distance of the Kp-Stations to the Auroral Zone
The average range of geomagnetic disturbance at subauroral latitudes shows a strong dependency on the distance to the auroral zone (e.g., figure 14 in Mayaud, 1980). The smaller the distance to the auroral zone, or the larger the CGM or QD latitude of a subauroral station, the larger is the average range of geomagnetic disturbance or K-variation. In K and Kp, this latitudinal dependency is compensated by a suitable choice of the K9-limits Mayaud, 1980). Geomagnetic secular variation and station replacements will change the QD latitude of the Kp-stations with time and with it the average range of geomagnetic disturbance. However, once a K9-limit is chosen for an observatory, this K9-limit usually is not changed any more. This could lead to a bias in the time series of the Kp-stations' K indices and in Kp, which we want to investigate for the period 1932-2020.  Figure 4a, the change in unsigned QD latitude with time due to secular variation and station replacement is shown color-coded for each of the 13 Kp-stations (using unsigned QD latitude makes sure that a poleward change in QD latitude has a positive sign for both the northern and southern hemisphere). The change in QD latitude was strongest (almost −5°) for OTT and FRD in North America.
The Kp-stations' weighted (according to their weight in Kp) mean change in unsigned QD latitude is plotted in Figure 4b and amounts to −1.7° from 1932 to 2020. The change due to the Kp-station replacements amounts to only −0.06°. Thus, stations replacement contributed only negligible and secular variation is the by far dominant process for the mean decrease in the unsigned QD latitude of the Kp-stations since 1932.
The auroral zone, i.e., the maximum of aurora sightings, is at around 67° geomagnetic latitude at least since the 19th century (Akasofu, 2003;Chapman & Bartels, 1940;Fritz, 1881) and if the QD latitude of the auroral zone stayed constant since 1932, then the mean decrease by 1.7° in the unsigned QD latitude would be a good estimate for the mean increase in auroral zone distance since 1932. Smith et al. (2017) found that the location of the auroral electrojet has changed from 1970 to 2014 in the same way as the QD latitudes have changed, confirming that the location of the auroral zone remained constant with respect to QD latitude. Changes in the dipole moment of the Earth's magnetic field could also affect the location of the auroral zone, but the observed change in dipole moment from 1932 to 2020, and its relationship to the location of the polar cap boundary (Cnossen et al., 2012), a proxy for the auroral zone location, yield a negligible southward movement of the auroral zone by 0.2°. Hence, the mean distance of the Kp-stations' auroral zone distance increased by about 1.7° QD latitude. To estimate the effect on the Kp time series of this change in auroral zone distance, we consider Thomsen (2004), who investigated the expansion of the auroral oval due to magnetospheric convection. Thomsen (2004), based on studies by Gussenhoven et al. (1981Gussenhoven et al. ( , 1983, found a linear relationship between the latitude of the equatorward auroral boundary (the ABI index) and Kp. Kp decreases by 0.4 units per 1° CGM latitude change of the auroral boundary toward the equator. The CGM latitudes in this range are almost identical to QD latitude and we will use QD latitude here. The relationship is linear and holds for 1 ≤ Kp ≤ 7 (see Figure 1 in Thomsen, 2004). However, the relationship describes the cumulative effect from the intensification of the current systems during auroral expansion on the one hand, and from the change in auroral boundary on the other hand (Thomsen, 2004). The relationship is thus an upper limit for the effect of auroral boundary movement, which is equivalent to decreasing the auroral zone distance to the MATZKA ET AL.  Kp-stations. This results in a (negative) bias in the Kp series due to a change in auroral zone distance since 1932 of not larger than 1.7°*0.4/° or about two thirds of a unit. In order to get a better estimate of the bias in Kp, we now estimate the effect of current intensification during auroral expansion on Kp to subtract it from this upper limit. Smith et al. (2017) investigated the global auroral ovals with magnetic data from satellites in near-Earth orbit. Their study is limited to the auroral electrojet currents for Kp ≤ 4 (but according to Table 3, this accounts for >90% of all Kp values) and yields higher QD latitudes for the auroral oval than the electron precipitation data investigated by Thomsen (2004). Smith et al. (2017) suggest (their figure 7) roughly a doubling of the auroral electrojet current strength for an increase in Kp by 2 for 0 ≤ Kp ≤ 4. For a doubling of the currents, the range of the K-variation would double and given the logarithmic scale for K ≤ 4, Kp would increase by about 1. Hence, about half of the increase in Kp comes from current intensification, and the other half is due to the approaching of the auroral oval during expansion. The results by Smith et al. (2017) allow separating the effect of auroral expansion and current intensification. They are about equally strong, but this remains a rough estimate since the auroral electrojet is only one of a number of current systems that affect K and Kp and while some current systems, like the field-aligned currents associated with the auroral electrojet, vary similarly to the auroral electrojet, others may vary differently. After subtracting this current intensification estimate from the upper limit by Thomsen (2004), our best estimate for the effect of changing auroral zone distance on Kp is 0.2 units per ° QD latitude and thus the bias in the Kp series since 1932 is in the order of minus one third of a unit, a number that is small compared to, e.g., the year-to-year variability of Kp shown in Figure 3a.
This negative bias in the Kp series since 1932 is due to the change in auroral oval distance due to secular variation. Other possible effects of secular variation on Kp, e.g., a change in geomagnetic declination, are not accounted for here and likely much smaller.

Nowcast Kp at GFZ
As described in section 2.2, the definitive Kp index is calculated from the K indices derived at each of the Kp-stations. Each Kp-station sends its K indices semimonthly to Niemegk, where the definitive Kp index is derived. In parallel, Niemegk also derives the so-called nowcast K indices and nowcast Kp index directly from near real-time magnetic data provided by the Kp-stations. GFZ distributes both the nowcast and the definitive Kp index.
At GFZ, the FMI algorithm (FMI stands for Finnish Meteorological Institute, where the algorithm was developed, see Sucksdorff et al., 1991) is used to determine the nowcast K values from 1 min data, which is provided by the institutes operating the 13 Kp-stations in near real-time (NRT). Data gaps of up to 15 min are linearly interpolated. The NRT data typically arrives with a delay of 5-10 min. Any incoming NRT data as well as any incoming delayed data (e.g., due to data transfer problems) triggers a recalculation of the affected nowcast K indices (and consequently the Kp index, see next paragraph). In these NRT operations, we calculate the first estimate of the K index in the current three-hour interval once 30 min of data (until mid of August 2020: 90 min of data) are available in this interval. At that time, the nowcast K reflects only geomagnetic activity that is included in the available data. Therefore, the range of the K-variation and consequently the nowcast K at this stage can be significantly smaller than the nowcast K calculated once all data of the three-hour interval are available (again, the same applies to Kp). For all previous intervals, K from a specific station is only given if at least 90 min of data are available for this interval. The FMI method estimates the quiet curve of any given UT day from the 1 min data of the previous day, the day itself, and the subsequent day. Therefore, nowcast K (and hence Kp) values can change with time until all data is available for these three days. To prepare the NRT data for use with the FMI method, we always fill in the missing future observatory data from the current and subsequent day with constant values equal to the latest 1 min data received. Within a three-hour interval, the nowcast K indices from the different Kp-stations are compared with each other. The largest K value is set invalid for the calculation of Kp if the difference between the largest and the second largest K value exceeds 3 units. This happens very seldom and up to now was caused by erroneous spikes or jumps in the 1 min NRT data resulting in erroneous K = 9 values during periods of significantly lower geomagnetic activity.
The first nowcast Kp index for an interval is published once the K values from at least 5 stations are available, which typically happens 35 min after the start of the three-hour interval. At that time, the nowcast Kp reflects only geomagnetic disturbance that is included in the available data and normally its value is significantly smaller than the value obtained when all data within the three-hour interval is available for calculating the individual K indices. The timeline for the calculation and distribution of Kp is analog to the timeline described for K in the paragraph above: With new data coming in, the Kp values get updated and once all data for the respective three-hour interval is available, Kp changes very little since only the determination of the quiet curve is affected by newly incoming data. The quiet curve is not changing any more once all data of the following day is available. However, late delivery of data from the Kp stations that fill in data gaps could retrospectively change our nowcast Kp values. Only with the monthly upload to the archive at GFZ Data Services (typically one week after the month in question, see section 4), the nowcast Kp values become fixed (even in the case of remaining data gaps). The calculation of nowcast Kp from nowcast K and Ks was slightly modified on August 1, 2020, and is described in the following paragraph.
We investigate differences between the nowcast and definitive Kp for 2019 by calculating the conditional and normalized occurrence of nowcast Kp given a certain definitive Kp. The normalization guarantees that the sum of all nowcast values for each definitive Kp amounts to 100%. For example, as shown in Figure 5a, for the condition of definitive Kp = 0o, there occurs a nowcast Kp = 0, 0+, and 1− in 73.2%, 26.1%, and 0.7% of the cases, respectively. Mostly, agreement between definitive and nowcast is around 70%. It is always above 55%, with the exception of Kp = 5−, where the nowcast Kp is underestimating the definitive Kp in 66% of the cases. Overall, the definitive and nowcast values differ mostly by one third of a unit, very rarely by two thirds (note the non-linear color scale in Figure 5a) and never by more than two thirds. A noteworthy discrepancy is seen for definitive Kp = 0o, where slightly more than 25% of nowcast Kp = 0+. There are two obvious reasons to explain the observed differences: First, the underlying NRT observatory data is preliminary, i.e., it can be incomplete (and some of our nowcast Kp products indicate the number of contributing observatories with available data), erroneous data might not be completely removed, and its calibration can be different from that of later data versions. Second, the FMI algorithm is used for all 13 stations when calculating nowcast K, whereas the observatory-dependent algorithms detailed in section 2.4.4 are used for calculating definitive K. To investigate the second effect, we calculate Kp from definitive geomagnetic observatory data from 1995 to 2017 using the FMI method for all stations and compare it to the definitive Kp. This FMI-based Kp frequency distribution (blue horizontal lines in Figure 5b) has 15% less Kp = 0 and 6% more Kp = 0+ than the definitive Kp (gray bars in Figure 5b). Also, it mostly overestimates the number of 0+ ≤ Kp < 3 and mostly underestimates the number of Kp ≥ 3o. This fits well to Figure 5c, where 31% of definitive Kp = 0 yield an FMI-based Kp = 0+. Also, around 3o ≤ Kp ≤ 6o, the FMI-based Kp tends to be one third of a unit smaller than definitive Kp in about 25%-30% of the cases (color code orange in Figure 5c) and one third of a unit larger in only about 10% of the cases (color code greenish).
As indicated by the similarity of Figures 5a and 5c, the larger part of the differences between definitive and nowcast Kp can be attributed to the FMI algorithm, and not to NRT data quality. To compensate this effect of the FMI algorithm on nowcast Kp, we slightly rescale the transformation from FMI-based Ks to Kp such that the distribution (red horizontal lines in Figure 5b) of the new rescaled FMI-based Kp follows more closely the distribution of the definitive Kp (gray bars). The details of rescaling are given in Table S4. Indeed, the rescaling improves the nowcast Kp, as confirmed by Figure 5d. In particular, only 23% of definitive Kp = 0 yield a rescaled FMI-based Kp = 0+ (as compared to 31% before rescaling). Also, the rescaled FMI-based Kp tends to be more equally distributed around the definitive Kp value in the region between 3o ≤ Kp ≤ 6o, with 15%-20% of the values (yellowish color in Figure 5c) being one third smaller or one third larger than the definitive Kp. This rescaled FMI-based Kp is produced and distributed as GFZ's nowcast Kp since August 1, 2020.
NOAA distributes a product resembling the nowcast Kp index and based on a subset of 8 of the 13 Kp-stations. The NOAA index is provided as a single-digit index, i.e., only full integer values are provided on a three-day plot (https://www.swpc.noaa.gov/products/planetary-k-index) and in files back to 1994 (ftp://ftp. swpc.noaa.gov/pub/indices/old_indices/). We also calculate a single-digit version of GFZ's nowcast Kp and compare it to the definitive Kp for 2019 in Figure 5e, while in Figure 5f, we compare the NOAA nowcast to the definitive Kp for 2019. Agreement with definitive Kp is usually better for the GFZ nowcast than for the NOAA nowcast, except for Kp = 5, where NOAA and GFZ nowcast correctly predicted 83% and 73% of the cases. As for the GFZ nowcast, the NOAA product usually overestimates Kp = 0, with predicting values of 1 in 43% of the cases (12% for GFZ, and GFZs nowcast is further improved from August 2020 due to rescaling).

Distribution, DOI, License and Acknowledgment of Kp
Here, we discuss the distribution of Kp. An overview of the data distribution channels is given in Table S5 for GFZ and in Table S6   The DOI for Kp and the derived indices is 10.5880/Kp.0001. The DOI is a persistent identifier and while the data set grows with time, no value published under this DOI will ever be overwritten. If necessity arises in the future to correct already published values, then the corrected data set will be published with a new DOI (e.g., 10.5880/Kp.0002). Older DOIs and data sets will still be available. For each DOI, an additional versioning mechanism will be available to document changes to the files such as format changes, which do not affect the integrity of the data. The data set comprises Kp, ap, Ap (nowcast since 2019 and definitive since 1932, in WDC-format yearly ASCII files) and the International Quiet and Disturbed Days (table-like yearly ASCII-files). These data sets are updated on a monthly basis, typically one week after the end of the month. The DOI links to GFZ Data Services through https://doi.org/10.5880/Kp.0001, and this is the original data set to which all other copies must be identical, including the ones described in the following paragraph. More details are found in Table S5.
Another, traditional, source of the Kp index is the Kp-website of GFZ (https://www.gfz-potsdam.de/en/ kp-index/). This page links to a number of data distribution services described in the following and in Table S5. The traditional Kp server (ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/) is updated half-monthly with definitive values and musical diagrams (also called "note-script", Bartels, 1949), and monthly and yearly with Kp frequency distributions and average ap values, respectively. In the WDC-format (format description on the FTP server), the Kp, ap, Ap, CP and C9 indices are accompanied by adjusted F10.7 values until 2006 and by sunspot numbers until 2014, when the new S N sunspot number was introduced (Clette & Lefèvre, 2016). As a new service, we provide now SN9C9-plots with the new S N sunspot number to replace the R9C9-plots. Sudden storm commencements (SSC, excluding sudden impulses, SI) are also listed and included in the musical diagrams. First, we use SSCs as detected at Niemegk observatory and, once available, replace these with the SSC list of the International Service on Rapid Magnetic Variations at Observatory Ebro (http://www.obsebre.es/en/rapid). The GFZ Kp-website also links to constantly updated digital and graphical nowcast values for the current and previous month, some indicating the number of stations available for the nowcast value, and a color-coded bar plot for the current and previous six days. A new service is an online generator and download function for this type of plots (as png, jpg, gif, eps, pdf) for dates back to January 1, 1932, with the most recent Kp values being nowcast values and earlier values being definitive values.
Another, new, source (ftp://ftp.gfz-potsdam.de/pub/home/obs/Kp_ap_Ap_SN_F107/) is meant to give experts and non-experts convenient, NRT capable access to a blend of the most current definitive and nowcast Kp indices and other relevant parameters. These files contain Kp, ap, Ap, the observed F10.7 (recommended for ionospheric and atmospheric studies), the adjusted F10.7 (recommended for solar studies), see Tapping (2013), and the international (daily total) sunspot number S N introduced in 2015 (Clette & Lefèvre, 2016). The blank-separated file format is easy to read in by code or as spreadsheets. Kp values are given as float with three decimals (e.g., 0.667 instead of 0.7 or 1−) and, next to the date, a linear time scale in days since January 01, 1932 is given as well as Bartels rotations. Both the starting time and the mean time of each three-hour interval is given, as Kp represents geomagnetic disturbance over the three-hour interval and the time stamp indicating the center of the interval might be more useful for certain applications. Two file versions exist. One version contains all data from one day in one line. The other version contains only the data from one three-hour interval per line and provides only Kp and ap. Yearly files and a file containing all data since 1932 are updated daily. A file for the most recent 30 days is updated continuously to support NRT operations. The Kp and sunspot number data is a combination of definitive data when available and nowcast data otherwise and regularly updated to replace nowcast with definitive values.
The Kp data set is published under the CC BY 4.0 license. Non-commercial and commercial use and redistribution are permitted under the conditions of the license, making the Kp universally available and usable. One of the conditions is attribution by citing the data publication (Matzka et al., 2021) and the present publication and referring to the data source GFZ. Table S6 lists a number of data redistributors that redistribute the Kp, including the WDCs for Geomagnetism in Copenhagen and in Kyoto and NASA.

Summary and Conclusions
The Kp index and derived products have proven their usefulness for over 70 years and are particularly important for space weather research and services. A DOI (Matzka et al., 2021) and the CC BY 4.0 license secure the data set's integrity and its future scientific as well as societal usability.
We review the production of the definitive Kp and how it changed in its history. The frequency distribution and fundamental properties of Kp are demonstrated on data from 1932 to 2019 (see Table 1 and 3, Figure 3). We stress the importance of the subauroral location of Kp-stations and investigate the change in distance between the Kp-stations and the auroral oval due to geomagnetic secular variation and station relocation. This distance is estimated to have increased by 1.7° QD latitude since 1932 (Figure 4b), almost all of it due to secular variation, while station relocation is negligible here. This results in a bias of the Kp series toward smaller Kp values with time. We estimate this bias, which developed from 1932 to 2020, to be smaller than two thirds of a unit and likely to be in the order of one third of a unit.
We describe the production of the nowcast Kp at GFZ for the first time. The first estimate is issued as early as 35 min (before August 2020 it was >90 min) after the start of the respective three-hour interval and it is updated continuously with new data coming in. Also, since mid-August 2020, the nowcast Kp index is rescaled to better match the definitive values. It yields the correct definitive Kp index in 60% of the cases and we expect that the remaining 40% are now equally distributed between nowcast values being one third of a unit larger (20% of cases) and one third of a unit smaller (20% of the cases) than the definitive Kp.

Data Availability Statement
Geomagnetic data was partly taken from INTERMAGNET. The total sunspot number data are provided by Royal Observatory of Belgium, Brussels and can be downloaded from the SILSO website (http://sidc.be/ silso/home). The solar wind data with one-minute time resolution were obtained from the NASA OMNI web database (http://omniweb.gsfc.nasa.gov/). The NOAA nowcast Kp values were taken from ftp://ftp. swpc.noaa.gov/pub/indices/old_indices/. The official Kp values used in this study are available at ftp:// datapub.gfz-potsdam.de/download/10.5880.Kp.0001. Calculation of QD was facilitated by https://github. com/aburrell/apexpy.