The potential to quantitatively reconstruct the interplanetary conditions of the previous 100 years from geomagnetic observations was foreseen by Russell (1975) and first attempted by Feynman and Crooker (1978) and Gringauz (1981). These early attempts used the aa geomagnetic index, initially constructed for 1868–1968 by Mayaud (1972) from observations made by geomagnetic observatories in southern England and in Australia. The aa index is the arithmetic mean of the northern and southern hemisphere subindices, aa _N and aa _S , both of which are constructed from the range (between minimum and maximum) of the irregular variation (i.e., after elimination of the regular daily variation), observed in 3-h intervals in either of the two horizontal components (northward or eastward, whichever gives the larger value). This range is then quantized, using quasi-logarithmic band limits that are specific to the observatory, into k indices, which are converted to a _k values using a standard scale. This procedure was first introduced by Bartels et al. (1939) and has been discussed in reviews by Lockwood et al. (2018a) and Lockwood et al. (2019a). In the northern hemisphere, the aa stations used are Greenwich (IAGA code name GRW), Abinger (ABN), and Hartland (HAD), which yield the k-index data sequences k _GRW , k _ABN , and k _HAD , respectively, and intercalibrated values are combined into aa _N . In the southern hemisphere, the sites are Melbourne (MEL), Toolangi (TOO), and Canberra (CNB), which yield the data sequences k _MEL , k _TOO , and k _CNB and correspondingly formed the southern hemisphere aa index, aa _S . After the work of Mayaud, the aa index has subsequently been continued to the present day.

The study by Feynman and Crooker (1978) used the derived correlation of aa with B s V S W 2 , where B _s is the southward component of the near-Earth interplanetary magnetic field (IMF) and V _SW is the solar wind velocity. Because B _s (and the IMF magnitude, B) was observed to be constant over solar cycle 20, Feynman and Crooker (1978) assumed that it remained constant and used aa to infer the long-term change in V _SW ; however, we now understand solar cycle 20 to be anomalous, and subsequent cycles have shown that the assumptions of constant B and/or B _s are not valid (Crooker and Gringauz, 1993). Stamper et al. (1999) demonstrated that on annual timescales, which are needed in these studies to average out seasonal ionospheric effects and other dipole tilt effects, the IMF orientation is an approximately constant factor, so the variation in B _s reflects that in the IMF strength B. Lockwood et al. (2017) and Lockwood (2022) have subsequently demonstrated that the assumption that the IMF orientation is constant is only accurate to within a 1σ error of 42% for 1-day means, but this error falls to 10.3% for 27-day means and to 4.9% for 1-year means. Hence, on annual timescales, the aa index depends on the IMF magnitude B and the solar wind speed V _SW to an accuracy of order 5%. The first reconstruction to separate the effects of B and V _SW on aa was carried out by Lockwood et al. (1999), who used the fact that the 27-day recurrence index from aa (Sargent, 1979) has a different (stronger) dependence on V _SW than aa itself. Lockwood et al. (1999) were aware that near-Earth interplanetary conditions are a local parameter, predominantly relating to the streamer belt, and so evaluated a global solar/heliospheric parameter, the open solar flux (OSF), F _S from B. This was achieved by using the Parker spiral theory (Parker, 1958) to derive the radial IMF component B _r from B and then utilizing the observation by the Ulysses spacecraft (Balogh et al., 1995; Smith and Balogh, 1995) that B _r in the heliosphere, averaged over IMF sectors of “toward” or “away” field polarity (B _r < 0 and B _r > 0, respectively), is independent of heliographic latitude. More specifically, Lockwood et al. (1999) used the physical explanation of that result by Suess and Smith (1996) and Suess et al. (1998), which predicts that |B _r | is independent of heliographic latitude. Note that Lockwood et al. reconstructed the signed OSF, F _S , the total away flux which, because of Maxwell’s equation ∇ ⋅ B ⃗ = 0 , equals the total toward flux, and the unsigned OSF is 2F _s . Here, we used signed OSF, and readers should be aware that some other studies used the unsigned OSF, which is larger by a factor of 2.

From the data available from Lockwood et al. (1999), it was evident that there had been a significant rise in the average values of F _S over the 20th century. The rise in OSF was reported to be a doubling, a factor that has been contested (e.g., Mursula et al., 2004; Svalgaard and Cliver, 2005). Next, various corrections have been required and implemented along with improvements to the method. The original reconstructions of Lockwood et al. (1999) found that the optimum rise in 11-year running, boxcar averages ⟨ F S ⟩ 11 (used to smooth out the solar cycles) was between 2.09 × 10¹⁴ Wb in 1901 and 5.42 × 10¹⁴ Wb in 1987. Hence, this is rise was by 159%. However, this initial estimate did not allow for the “excess flux” (caused by the effects on |B _r | of “folds” or “switchbacks” in the heliospheric field and by the orthogardenhose flux; see discussion below). Lockwood et al. (2014b) used the Monte Carlo fitting technique and four pairings of indices to generate an ensemble of estimates that made allowance for the excess flux. These reconstructions gave a minimum ⟨ F S ⟩ 11 of 2.06 × 10¹⁴ Wb in 1901 and an earlier and smaller maximum of 4.51 × 10¹⁴ Wb in 1955, and hence an estimated rise of 110% with a 2σ uncertainty range between 95% and 125%. The uncertainty range is large because of the uncertainty in the minimum values in solar cycle 14 at the start of the 20th century: very small changes to this value make disproportionately large changes to the percentage rise estimate.

Lockwood (2001), Lockwood (2003), and Lockwood et al. (2009b) showed that the OSF subsequently had reached a peak and was falling again and, as discussed later in the present study, both geomagnetic and direct in situ spacecraft observations show that ⟨ F S ⟩ 11 fell back to 2.64 × 10¹⁵ Wb by 2009, at the start of solar cycle 24. This was the lowest value since 1902, at the start of solar cycle 14. Furthermore, computation of the heliospheric modulation potential and cosmic ray fluxes from the F _S variation revealed that the twin peaks formed a grand solar maximum (Usoskin et al., 2021). Steinhilber et al. (2012) used cosmogenic isotopes to estimate the near-Earth IMF B (which is related to OSF but is only an approximate proxy for OSF because of variations of the Parker spiral associated with solar wind speed) and found that the peak value in the 20th century was of a magnitude that has been seen just 24 times previously in the 9000-year cosmogenic isotope record, an average repeat period (of peaks in the IMF magnitude equal to or exceeding the 20th-century peak) of 360 years (Barnard et al., 2011). Indeed, we now know that the interval since 1610, when Thomas Herriot made the first telescopic observations of sunspots, contains a grand solar minimum, the Maunder minimum (Usoskin et al., 2015), and a grand solar maximum (Lockwood et al., 2009b) and covers almost the full range of solar activity levels seen over the approximately 12,000-years interval over which cosmogenic isotope data can give information about solar activity (Usoskin, 2017). At the start of solar cycle 24, Barnard et al. (2011) investigated the rise and fall of the 20th century grand maximum and compared it with 24 previous analogous maxima deduced from cosmogenic isotope abundance measurements for the past 9000 years. They made analogue forecasts that predicted a probability of approximately 8% that the fall would be rapid enough for grand minimum conditions to commence around 2050. However, the optimum prediction (the median of the ensemble of 24) was for a continuing slower decline than this (in particular, at about the same rate as seen over cycles 22 and 23) in sunspot numbers and OSF, which is broadly consistent with what has been subsequently observed in solar cycle 24.

Because this 20th-century grand maximum peaked in the interval of modern spacecraft data, we here refer to it as the “modern grand solar maximum,” hereafter “MGSM.” There is great interest in this grand maximum because of its potential to tell us about solar evolution, long-term variations in space weather, and the solar magnetic cycle (see studies in the collection edited by Banerjee et al., 2018).

The MGSM can be seen in several datasets that extend back far enough in time and that continue to the present day. Figure 1 shows two examples, sunspot numbers and the homogeneous aa index, aa _H . The top panel shows the newly recalibrated international sunspot number, R _ISN (Clette and Lefèvre, 2016), and the revised sunspot group number, R _G (Vaquero et al., 2016). The latter is scaled by a factor f, the ratio of the means over the interval 1978–2018.

The aa _H index is based on the same observations as Mayaud’s classic aa index, but the long-term drifts in the magnetic latitudes of the stations are allowed for by correcting for the distances between the stations and an average auroral oval location. This is done using the IGRF model and (for before 1900) the gufm1 model of the intrinsic geomagnetic field. This is a considerable advance on the original aa index for which these distances were used but assumed constant for each station (Lockwood et al., 2018a). In addition, improved intercalibrations (that allow for the time of year) between data from different stations were used, along with a model to remove the time-of-year/time-of-day effects at a given station location (Lockwood et al., 2018b). These changes result in aa _H , unlike aa, having two hemispheric subindices (aa _HN and aa _HS for the northern and southern hemispheres, respectively) with almost identical long-term variations, despite them being entirely independent in their derivation. In other words, the much better agreement of aa _HN and aa _HS (compared with that for aa _N and aa _S of the original aa index) was achieved from the improvements to the derivation algorithm and not by any intercalibration to bring them closer together. The aa _H index data series runs from 1868 to the present day and here has been extended back to 1844 using the k-index data from the Helsinki and St. Petersburg magnetometers, using the procedure described by Lockwood et al. (2013a) and Lockwood et al. (2014b). The history of the development of this extension is reviewed in Section 4 of Lockwood et al. (2014b).

In both panels of Figure 1, the thinner lines are annual means, and the thicker lines are 11-point boxcar (running) means of the annual data. After about 1925, the behaviors of R _G and R _ISN are very similar. Sunspot numbers have minima that return to close to zero in every cycle, although there is some cycle-to-cycle variation in the minima as well as in the maxima. The 11-year running means show two clear peaks, the first being the larger of the two because of the exceptionally large amplitude of solar cycle 19 (June 1954 to October 1964). The 11-year averages of cycle 24 (December 2008 to December 2018) are almost exactly the same as they had been in cycle 14 (January 1902 to July 1913) with ⟨ R ISN ⟩ 11yr ≈ 50. The peak value in cycle 19 was more than double this value ( ⟨ R ISN ⟩ 11yr ≈ 125), so both the rise and fall in average sunspot number over the period of 1900–2020 are by 150%. In view of this, a large rise and fall in OSF is not surprising. The behavior of the aa _H index also shows a double peak, but differences exist. It is well known that the minima in annual aa and aa _H values show much greater cycle-to-cycle change than the corresponding minima for sunspot numbers. In addition, the peaks in ⟨ a a H ⟩ 11yr are more triangular in shape than for the sunspot numbers, and the second peak is slightly larger than the two. The two peaks in the 11-year running means between weak solar cycles 14 and 24 define what we term here as the MGSM.

Figure 2 shows this grand maximum in two other long solar–terrestrial observation sequences in the same format as Figure 1, both of which commence in the early 1930s. The upper panel is the ionospheric F2 peak critical frequency, foF2, from the intercalibrated Slough and Chilton ionosondes. Near-noon values are shown, being the means of values taken at 11, 12, and 13 UT. The double peak is similar to that in sunspot numbers, and the variation is consistent with that in R _ISN , given the nonlinear regression of annual means of these foF2 data and R _ISN , as demonstrated by Lockwood et al. (2016b). Note that the fractional change in foF2 is smaller than that for R _ISN because the annual means of foF2 fall to minima of approximately 5 MHz, whereas R _ISN falls to minima near zero. The bottom panel is for the planetary ap geomagnetic range index. This is quite similar to the plot for aa _H given in Figure 1A. However, there are some important differences. In particular, the first of the two peaks of the MGSM is the larger for ⟨ap⟩₁₁, whereas for ⟨ a a H ⟩ 11 , the second peak is larger. Because ap is based on a planetary network of stations and aa _H is based on just two stations (currently Hartland (HAD) and Canberra (CNB)), this initially appears to point to a problem with the calibration of the aa _H index. However, in this study, we show that aa _H is in excellent agreement with the k-index data from the Niemegk (NGK) station, k _NGK . This is highly significant because ap is constructed by rescaling the k indices from the other stations in the ap network to what Niemegk would have seen using standard lookup tables before averaging (Lockwood et al., 2019a). Hence, this difference between aa _H and ap appears to be arising from the procedure used to compile ap and from the fact there have been numerous changes to the network of stations used to compile the ap index (Lockwood et al., 2019a). This stresses a key point that the homogeneity in the construction of a geomagnetic index is vital if we wish to use it to reconstruct the space weather conditions through the MGSM.

It is self-evident that the reconstruction of past solar and interplanetary conditions from historic geomagnetic data critically depends on the stability of the calibration of those data. There have been a number of instances of calibration skips and drifts in data from a given station that have emerged from comparisons with data from different stations. For example, in developing the IHV index (Svalgaard et al., 2004; Svalgaard and Cliver, 2007) using data from many stations, Leif Svalgaard in 2004 noted in a letter to the data curators that there appeared to be an error in the early hourly mean data from the Eskdalemuir station, something also noted by Martini and Mursula (2006). From careful comparisons with data from nearby stations, this error was corrected by Macmillan and Clarke (2011). Another example was poor calibration of the horizontal force variometer at Helsinki for a 6-year interval noted by Svalgaard (2014), which was causing the horizontal H component to be too small; Lockwood et al. (2014a) also found that correcting for this using data on the vertical Z component indeed brought the data in line with the observations from other stations, particularly at the nearby St. Petersburg station.

Such corrections have been very valuable and important. However, there is a potential and important pitfall here of which we need to remain aware or we will make erroneous corrections that would be a serious retrograde step. The key point is that if the data sequences from two stations diverge, there are a number of possibilities: first, one of the data sequences is in error (in which case, which one?); second, they could both be in error; third, they could both be correct, and it is the expectation that they should be the same that is in error. The first and second possibilities are self-evident, but the third is more subtle. We know that the geomagnetic activity response depends on the geomagnetic latitude of the station. The secular changes in Earth’s intrinsic magnetic field mean that data from different stations will drift relative to one another (see Lockwood et al., 2018a), and we must be careful not to falsely ascribe these changes to measurement error. The reason is that making a correction based on the assumption that the difference is because of an instrumental error effect means that the variation will be in error when we apply a correction for the secular drift of the station. Furthermore, the sensitivities of the geomagnetic observations to variations in the IMF, B, and solar wind speed, V _SW , depend on the station location. Therefore, correcting data from two stations to make them agree, on the false assumption that they should agree, will cause errors in both the derived B and V _SW . Another potential pitfall of making corrections is the wrong assumption that the geomagnetic activity data from different stations are linearly proportional because, in reality, local site conditions may mean this is not the case.

As expected, if we have one station that disagrees with several others in similar locations on a given measurement, it is an indication that it has a calibration problem, and this is best dealt with by using the median of the distribution of all the available values (with it being less prone to distortion by poorly calibrated stations than the mean value). The problem is that as we go back in time, we have fewer operating stations, and so the uncertainties associated with using the ensemble of all available stations inevitably increase. Furthermore, the geographic distribution of those stations changes, which introduces an inhomogeneity into the data series. For this reason, the second approach is to use only a small number of stations but try to make the data from them as homogeneous as possible (including allowing for the secular drift in geomagnetic latitude and site differences). Hence, there are two philosophies to generating a reliable geomagnetic data sequence to extend modern spacecraft data on the interplanetary medium back in time: (A) take the mean, or preferably the median, of the distribution of all the stations available at a given time and (B) use a limited number of stations but do all one can to ensure the data series is homogeneous in the sense that data from the past is consistent with the data from the space age (the latter being the data that are compared with interplanetary observations). Approach A has the advantage of using all the available data but the disadvantage that the number of stations decreases as we go back in time, thereby increasing the uncertainties. This also means that the data for early years is not from the same network of station locations as that during the space age. When the number of stations is high, such an index is genuinely global, but as one goes back in time, it becomes increasingly local, applicable only to where the few early stations were sited and not global in nature. Hence, it is not true that using all the data available is necessarily the best approach. Approach B could, in theory, achieve uncertainties that are close to being constant throughout the data series and at all times using data obtained from stations reasonably close to the stations giving the data during the space age. Philosophy B uses model values of the intrinsic magnetic field to allow for the (changing) magnetic latitude of stations and so turn the local measurements into a global estimate. However, it is not deploying all the data available and so is more prone to errors introduced by station calibration errors.

The IDV (interdiurnal variability) geomagnetic index (Svalgaard and Cliver, 2005) and IHV (interhour variability) geomagnetic index (Svalgaard et al., 2004; Svalgaard and Cliver, 2007) are examples of the application of philosophy A, whereas the aa _H (Lockwood et al., 2018a,b) and IDV(1d) (Lockwood et al., 2013a,b, 2014a) indices are examples of the application of philosophy B.

Given that geomagnetic responses vary with location, we should not attempt to make all the geomagnetic data agree, because that is based on the incorrect assumption that they should all agree. However, when we use geomagnetic data to reconstruct interplanetary conditions, the results should not depend on where on the surface of the Earth the data come from; hence comparing the interplanetary reconstructions obtained using different geomagnetic data sequences is the most valid test.

Recent studies have confirmed that there is no such thing as a “universal coupling function,” a combination of interplanetary parameters that accurately predicts all terrestrial responses to interplanetary variations. We know this because even different geomagnetic activity indices are found to respond differently to a given set of interplanetary conditions, let alone other space weather activity indicators (Lockwood and McWilliams, 2021; Lockwood, 2022).

Table 1 presents the results of a multiple regression analysis of the four indices used in the present study, aa _H , IHV, IDV(1d), and IDV. We used the MATLAB two-dimensional polynomial fitting routine “fit2dPolySVD” and evaluated the polynomial using “eval2dPoly” (Whitehead, 2021). We use polynomials of order 2 and, based on previous studies, fit using the near-Earth IMF B and solar wind speed V _SW . Table 1 presents the six coefficients (a to f) that yield the predicted index, I _p1, given by I p 1 = a + b . V S W + c . V S W 2 + d . B + e . B . V S W + f . B 2 (1)

We then evaluate the fit by computing the correlation coefficient between the sequence of I _p1 values from the annual means of B and V _SW from the Omni interplanetary data (King and Papitashvili, 2005) and the annual means of the observed index I _obs . The optimum correlation coefficient from this procedure is r _opt1 , and the 2σ uncertainty range in that value is between r _min1 and r _max1 .

The left-hand side of Equation 1 has just two free variables, B and V _SW , and so in theory, we could solve for both using just two geomagnetic indices. However, the combination of terms in Equation 1 means the solution is complex and almost certainly would require a numerical iteration technique. Table 1 demonstrates that we obtain almost as high correlations with the much simpler formulation: I p 2 = g + h . B . V S W n (2) which was used by Lockwood et al. (2014b) and is much more readily solved. As for Equation 1, we evaluate the minimum, optimum, and maximum correlation coefficients between I _p2 and I _obs , which were r _min2 , r _opt2 , and r _max2 , respectively. All these correlation coefficients are presented in Table 1 for each of the four I _obs indices used, along with the 2σ errors of the optimum correlations in both cases (ϵ ₁ and ϵ ₂) and the differences between the correlations for the two fit procedures Δ _min = (r _min2 − r _min1 ), Δ _opt = (r _opt2 − r _opt1 ), and Δ _max = (r _max2 − r _max1 ). Table 1 shows that all the correlations are exceptionally high, with the lowest optimum value being 0.915 and the highest being 0.971, but that the Δ values are positive, showing that using Equation 1 always generates better fits than using Equation 2, as we would expect, because it has six free fit parameters, as opposed to the three in Equation 2. However, these Δ values are very small and considerably smaller than the 2σ uncertainties, ϵ ₁ and ϵ ₂. We also applied the Meng-Z test for the significance between two correlations (allowing for the intercorrelation of the parameters) (Meng et al., 1992) and p-values for the null hypothesis that there is no significant difference between the two correlations were always above 0.85. Furthermore, we used fit residual analysis to determine the AIC of each of the fits, and again, in no case did the fit using Equation 1 give a significantly lower value than the fit using Equation 2. These tests justify using the much simpler fit procedure with fewer free fit parameters (Equation. 2).

The procedure used is a variant of that used by Lockwood et al. (2014b), but with an additional solar cycle of data now available (except in the case of IHV, which has not been updated). That extra cycle (solar cycle number 24) is much weaker than the others previously available in the space age and so extends the range of values covered by the fits, thereby reducing the degree of extrapolation involved in the reconstruction of the lower values derived for around the start of the 20th century. Figure 6 shows the variations of B . V S W n (for B in nT and V _SW in kms ⁻¹) for the optimum exponent n and the best-fit annual means of the observed geomagnetic index, I _obs , from Equation 2, (I _obs − g)/h, for the four geomagnetic indices discussed and checked in Section 3.

The reconstruction procedure is the most sophisticated yet used as it is the only one (as far as we are aware) that employs a Monte Carlo technique, based on the observed fit uncertainties, to generate an ensemble of 1000 reconstructions of B . V S W n for a given index. The root-mean-square deviations of the best fits and the data ϵ o 2 = ⟨ ( I obs − I p 2 ) 2 ⟩ 1 / 2 are then computed. To generate each ensemble member, a random number generator was used to perturb all the data points individually by an error drawn at random from a Gaussian distribution of mean zero and standard deviation ϵ _o2. The best fit is then found for each set of error-perturbed data using the Nelder–Mead search method, and we check that 95% of the unperturbed data points lie within the plus and the minus 2σ points of the distribution of 1000 fitted variations. We then use Equations 4–7 of Lockwood et al. (2014b) but in a slightly different way in that we generate an ensemble of 1 million members for each of the four index pairings by using every permutation of the 1000-fit ensembles for the pair of indices used. Because of the matrix manipulation used in MATLAB, this can be achieved in relatively short computation times. This is done for each of the four pairings of indices that show sufficiently different best-fit values of the exponent n. The annual median of the total of four million reconstruction ensemble members is then taken for both B and V _SW , along with the 2σ uncertainty limits.

Figure 7 studies how well the exponent n is defined in each case. The top panel, A, shows the correlation coefficients r ₂ between the index I _obs and B . V S W n as a function of the exponent n. The n values giving peak correlation r ₂ = r _opt2 are marked by vertical dashed lines. The question arises as to what extent small differences in r ₂ are statistically significant. We here look at the significance S of the difference between r _opt2 and r ₂ at a different n by computing the p-value of the null hypothesis that they are the same (where S = 1 − p). To compute p, we use the Meng-Z test for the difference between the correlations between A and B and between A and C, which allows for the intercorrelation of B and C (Meng et al., 1992). We test against the AR1 red-noise model by using the effective number of independent data pairs, N _eff , given by N eff = N 1 − a 1 / 1 + a 1 (3) where N is the actual number of data pairs, and a ₁ is the autocorrelation at lag 1 (Wilks, 1995). The results for S as a function of n are shown in the middle panel (B) of Figure 7.

If two indices have identical exponents n, there is no independent information in the two, and B and V _SW cannot be separated—the errors in their derivation become infinite. If the n values are similar such that the difference approaches zero, the errors in B and V _SW are very large and tend to infinity. For the four indices used here, there are six permutations for combining them in pairs. For two index pairs, specifically IDV and IDV(1d) and IHV and aa _H , the peak correlations are at almost identical n, which means they cannot be used to separate the variations of B and V _SW . This leaves four usable permutations. The bottom panel (C) of Figure 7 checks that the values of n for these indices are significantly different. The probability of index I _obs1 having a general value of n is (1 − S ₁), and the probability that a second index I _obs2 has the same value is (1 − S ₂); hence, the probability that both have that value of n is (1 − S ₁) (1 − S ₂), which is plotted as a function of n in the bottom panel. It can be seen that the most clearly distinct combination is IDV and aa _H (in red) for which the peak probability of them not being distinct is 0.004% for n near 1. The least-distinct pairing is IHV and IDV(1d) (in orange) for which the peak probability is larger (0.08%) but still very small. For IHV and IDV (in green), it is 0.07%, and for IDV(1d) and aa _H (in black), it is 0.06%. The uncertainties in B and V _SW from a given index pairing are found to increase with the peak value of (1 − S ₁)×(1 − S ₂) and are extremely large, and indeed often infinite, for the two unused pairings, which give (1 − S ₁)×(1 − S ₂) that peaks very close to unity around the peaks of the correlations.

Figure 8A shows the reconstructions of the near-Earth IMF, B. The four colored lines are the median values of the ensembles of 1 million fits for each of the four index pairings. The blue dots give the annual means of the observations by spacecraft in near-Earth interplanetary space from the Omni dataset, B′. Figure 8B presents the superposed scatter plots of B as a function of B′ for 1964–2020, inclusive, and Table 2 gives the correlation coefficients. The four best-fit least-squares linear regression fits are shown in Figure 8B and are indistinguishable. For all years after 1880, the four median reconstructions are almost identical and lying very close to the center of the gray area, which is the ±2σ uncertainty band of the ensemble of all four million values. This is the ultimate test of the geomagnetic data used to compile all four reconstructions because the interplanetary parameters do not depend on where the geomagnetic data were measured, nor on the instrumentation used, nor on the calibration of the instrumentation. Hence, if there are differences in the geomagnetic data that are used, they should be there and should not be corrected on the (false) assumption that the geomagnetic data should all agree with each other.

Before 1880, the situation is different. In this interval, the only data available are the aa _H index (extended back to before 1868 using data from Helsinki, St. Petersburg, and, where possible, Greenwich), IDV(1d) and IDV (Lockwood et al., 2014a). As discussed above, IDV(1d) at this time is based on the Helsinki data (after correction for erroneous hourly means of the horizontal component during cycle 11), whereas IDV uses Bartels’ u-index, which he regarded as illustrative only for most of the interval. The aa _H –IDV(1d) pairing gives the black line and the aa _H –IDV pairing the red line, and we can see that before 1880, the two give quite similar variations but do diverge somewhat. The procedure at this time generates ensembles of 2 million members: the median of the total ensemble is between the red and black lines, in which the ±2σ uncertainty band is somewhat broader than that for later years.

Figure 9 shows the corresponding four median reconstructions of the near-Earth solar wind speed, V _SW . Again, the annual mean from the Omni dataset of spacecraft observations V S W ′ is given by the blue dots; the scatter plots of V _SW as a function of V S W ′ are given in Part B and the correlation coefficients in Table 2. The agreement between the different pairings is not as close as for the IMF B, nor with the spacecraft observations for 1964–2020. However, agreement is still very good. All the time series show peaks in annual means of V _SW in the declining phase of most solar (but not all) cycles. This is expected because in this phase of the solar cycle, isolated coronal holes form at low latitudes, along with extensions of polar coronal holes to low latitudes. This increases the probability of Earth-directed fast solar wind and raises the occurrence of corotating interaction regions impacting Earth, giving recurrent geomagnetic disturbances. The derived sequence has similarities to that derived by Mursula et al. (2017) using the fact that the geomagnetic activity response depends on the latitude of the station (Finch et al., 2008). Mursula et al. (2017) noted variations in this time series related to the rise and fall of the grand maximum. Here, we again note that the close agreement between the results of different pairings of geomagnetic activity indices means that errors in those indices are small, and so, differences are real and should not be corrected.

We here use a development of the procedure of Lockwood et al. (2014b) to compute the (signed) OSF F _S . We use the Parker spiral theory, and we need to evaluate its accuracy and make any correction required.

Figure 10 shows the results of a survey of daily mean values of the interplanetary data and compares the observed IMF “gardenhose” angle with those predicted by the Parker spiral theory. Polar histograms of the observed (gray filled, θ _obs ) and predicted (mauve outlined, θ _p ) distributions of the IMF gardenhose angle are presented for daily means of the near-continuous interplanetary observations from 1995 to 2020, where θ obs = t a n − 1 − B Y / B X (4) and B _Y and B _X are the Y and X IMF components in the GSE frame, respectively. The plots are for 12 ranges of observed solar wind speed, V _SW , each containing one 12th of the total number of samples in the dataset. The predictions use the daily mean V _SW value and Parker spiral theory θ p = t a n − 1 r ω / V S W (5) where ω is the angular velocity of the solar corona with respect to the fixed stars, and r is heliocentric distance. A corresponding plot for hourly means was presented by Lockwood et al. (2019), who studied the causes of the largest deviations from Parker spiral orientations, namely, orthogardenhose flux. We predict the modulus of the radial field | [ B r ] p | using ∣ B r p ∣ = B ⁡ c o s β ≈ B ⁡ c o s θ p = B r ω / r 2 ω 2 + V S W 2 1 / 2 (6) where β is the angle between the IMF vector B ⃗ and the radial direction. This is, in general, different from the gardenhose angle θ, and the two are only the same if the IMF component B _Z = 0. Hence, we need to evaluate the effect of nonzero B _Z caused by the need to adopt the approximation inherent in Equation 5 because we cannot reconstruct B _Z . We also need to investigate the effect of the averaging timescale τ as we will be using annual mean data. The combined effect of these two factors is investigated in Figure 11, compiled from the 26 years’ data used to compile Figure 10. We evaluate ∣ [ B r ] p ∣ using Equation 6 and from that the fractional error given by ϵ = | B r p | − | B r obs | / | B r obs | = B ⁡ c o s θ p / | B r obs | − 1 (7) where | [ B r ] obs | is the absolute value of the observed radial IMF component. The distribution of ϵ values was computed for 1-day means (shown in Figure 11A) and then for running boxcar means of the 1-day values over timescales τ of B 7 days, C 27 days, and D 1 year. Figure 11D shows that the mean of the distribution of errors μ _ϵ is not zero, as was assumed by Lockwood et al. (2014b), but is 0.198, which means that we need to divide the predicted radial field estimates | [ B r ] p | by (1 + μ _ϵ ) ≈ 1.2. We find no evidence that the error ϵ depends on either | [ B r ] p | or |[Br] _obs |, and so the same correction can be applied to all data. This correction explains why F _S estimates are a bit smaller than those given by Lockwood et al. (2014b) but do not change the estimated percentage rise because all F _S values are changed by the same factor.

The mean μ _ϵ and standard deviation σ _ϵ of the distribution of ϵ values are given in each panel along with the 1σ fractional error in using the procedure and then correcting for the offset of the mean value, σ _ϵ /(1 + μ _ϵ ). For daily means, that error is by a factor of close to 6, but the distribution narrows as τ is increased, such that for τ = 7 days, the error is 25%; for τ = 27 days, it is 12.5%; and for τ = 1 year, it is just below 5%.

We also have to consider another factor that has been given a variety of names: “excess flux” (Lockwood et al., 2009b,a; Lockwood and Owens, 2009), “inversions” (Owens et al., 2013), “folded flux” (Owens et al., 2017b), and “switchbacks” (e.g., Mozer et al., 2020). The radial heliospheric magnetic field (HMF) can be used to compute the OSF because of two factors. First, we are averaging out any structure in the Carrington longitude by taking means over several solar rotations for annual means at Earth (specifically 13.42 rotations in a leap year and 13.38 in a non-leap year). The second factor is that the Ulysses satellite showed that the radial component of the heliospheric field was approximately independent of heliographic latitude (Balogh et al., 1995; Smith and Balogh, 1995), a result we can understand considering that close to the Sun, the plasma beta is low, i.e., pressure is dominated by the magnetic field. This means solar wind flows close to the Sun are slightly nonradial until the tangential pressure is equalized, which means that the modulus of the radial magnetic field |B _r | is equalized, and this is carried out by the radial solar wind flow into the heliosphere (Suess and Smith, 1996; Suess et al., 1996, 1998). The tangential magnetic pressure depends on the radial field of either polarity (toward or away from the Sun), and this is why it is the modulus of the radial field that is constant. At a heliocentric distance r, the unsigned magnetic flux (of either polarity) threading the surface of radius r is 4πr ²|B _r |. Owens et al. (2008) surveyed radial field measurements throughout the heliosphere and showed that although |B _r | remained approximately independent of latitude at all r, the unsigned flux F _S increased with r. We here use signed flux (of one polarity), and Maxwell’s equation ( ∇ . B ⃗ = 0 (the nonexistence of magnetic monopoles) means that for any r, the outward magnetic flux equals the inward flux, and the signed flux is F _S = 2πr ²|B _r |. At the coronal source surface r = r _o , where B _r = B _ro , this flux is the signed OSF F S = 2 π r o 2 | B r o | , and from the result of Owens et al. (2008), we can write F S = 2 π r o 2 | B r o | = 2 π r 2 | B r | − Δ r (8) where Δ is the “excess flux” that increases with r (Lockwood and Owens, 2009). This is the magnetic flux that threads the coronal source surface (with either “Toward,” T, or “Away,” A, field polarity in the − r and + r directions, respectively) and loops back toward the Sun in the heliosphere (termed an “inversion,” a “switchback,” or “folded flux”) so that it threads a surface at a given heliocentric distance r more than once. Excess flux also includes orthogardenhose flux tubes that are also bent back toward the Sun but by a smaller angle, such that they are no longer in the gardenhose sectors of IMF orientation. If it does this once, such a heliospheric flux tube containing OSF flux dF _S of A polarity at the coronal source surface contributes flux 2. dF _S to total A flux threading the sphere of radius r and adds dF _S to the total T flux threading the sphere. Likewise, if the flux dF _S at the source surface is of T polarity, it contributes 2. dF _S to the total T flux threading the sphere of radius r and dF _S to the total A flux threading the sphere. In both cases, the excess flux added is Δ = 2. dF _S (Lockwood and Owens, 2013). If the flux tube were to do this twice, it would add Δ = 4. dF _S , and so on. Note that excess flux is undoubtedly a real physical phenomenon and cannot be dismissed as noise or an artifact of using a modulus—indeed, the use of the modulus is necessary here to avoid artifacts induced by averaging over a region containing two opposite polarities. Note that the theory of why the heliospheric radial field is constant (Suess and Smith, 1996; Suess et al., 1996, 1998) shows that it is the modulus that we should use because both the T and A fields contribute to the tangential field pressure. Note also that excess flux has two causes. The first is orthogardenhose flux, where flux tubes are bent back toward the Sun but by an angle large enough to take it out of the gardenhose IMF orientation sector but not large enough to take it into the gardenhose sector of the opposite T/A polarity. The second is folded flux (a.k.a. switchbacks), where the flux tube is bent back by an angle large enough (approaching 180°) for it to remain in a gardenhose orientation sector (but becomes T flux where it was A or vice versa). Both contribute to excess flux, but only in the orthogardenhose flux case can it be identified from the field orientation. Recent results from the Parker Solar Probe mission showed how common switchbacks are in the inner heliosphere (e.g., Mozer et al., 2020), and so accurately computing the excess flux Δ is of great importance.

However, there is no easy way to compute the total excess flux, Δ. From potential field source surface (PFSS) modeling, we can identify the polarity of the radial field at the source surface and average over the T and A sectors of the source surface. If we could identify where sector boundaries in the source surface map to in the satellite data, we could average the signed B _r over the intervals of one source polarity (it is important to note that we need source polarity, not polarity at the spacecraft), which would cancel excess flux because the additional + dF _S of A flux of any orthogardenhose or folded flux tube would cancel the −dF _S of the corresponding T flux. Smith and Balogh (1995) first identified the latitudinal invariance of the radial field (often referred to as the “Ulysses result”) by assuming that they could define the source sector boundaries from the field polarity at the spacecraft. However, the T/A polarity at the source surface does not always stay the same all the way to the spacecraft, and so, we cannot, in general, identify the source sector boundaries in the satellite data. For example, if there is a short interval of opposite polarity radial field in a sector, one can never know if it is folded flux (and so does not appear in the source surface) or if it reflects real polarity structure on the source surface. Hence, this method does not provide a self-consistent method, and we cannot rely on this method to average out Δ (Lockwood and Owens, 2009). In practice, comparison of OSF values from the (PFSS) modeling of the corona can be made to agree by simply averaging the radial heliospheric field at the spacecraft over a given timescale. Some of these agreements are achieved by averaging out toward and away excess flux, but toward and away structure that maps back to the source surface, and hence is genuine OSF, is also averaged out. It was found that an averaging timescale of approximately 1 day made spacecraft and PFSS estimates agree rather well, but there was no physical reason to choose this averaging timescale, and so, it was an arbitrary correction, and we could never be sure it always applied. In addition, it must be remembered that the PFSS method contains many assumptions.

Other methods, such as using orientation with respect to the Parker spiral theory (Erdős and Balogh, 2012, 2014), can be used to remove the orthogardenhose component of the excess flux but relied on the use of an averaging timescale over an arbitrary timescale to remove the folded flux component and so were not any more satisfactory. Lockwood et al. (2009a) devised a method to use the gradients in solar wind speed with frozen-in flux to estimate the excess flux (the kinematic correction) and showed that it could broadly explain the difference between PFSS and spacecraft values more physically. However, by far, the best method is that devised by Owens et al. (2017b) using strahl electrons to identify folded flux (and orthogardenhose flux) and compute the excess flux. Strahl electrons are generated by the high temperatures of the solar corona. They flow along heliospheric field lines, and because the field-aligned velocity is much higher than field-perpendicular velocities of the field lines, the field lines move only a small distance during the short electron transit times, and the electron trajectories are close to (but not quite exactly) field-aligned. Strahl electrons are sometimes seen moving back toward the Sun, and this is a unique identifier of excess flux (either orthogardenhose or folded flux). Strahl on orthogardenhose HMF is sunward, as is strahl on inverted gardenhose HMF that has been through a switchback and is folded back on itself. Owens et al. (2017b) used strahl data from the Advanced Composition Explorer (ACE) spacecraft to compute the total flux of field lines showing sunward strahl and properly measure the excess flux Δ for the first time. The results showed that the kinematic correction was valuable but slightly and systematically overestimated the excess flux. Frost et al. (2022) have recently extended the survey of Owens et al. (2017b) by using both the ACE and WIND spacecraft, and we here use linear regression their values of the OSF derived using Equation 8 to scale OSF values obtained from the geomagnetic reconstructions of B and V _SW using Equation 6. This yields F S ′ = s F S G + c = s 2 π r 2 | B r p | + c = 2 s π r 2 B V sω / r 2 ω 2 + V S W 2 1 / 2 + c (9) where s and c are the coefficients from the linear regression of the estimates F S ′ from spacecraft data (allowing for excess flux using the strahl observations), and F _SG is the estimate derived from the annual mean geomagnetic activity (with no allowance for excess flux). Note that the excess flux Δ is accounted for by the intercept of the regression, c and the slope s allows for the mean fractional error μ _ϵ introduced by using the Parker spiral theory for annual means (see Figure 11). Table 2 presents the coefficients s and c for the OSF values F S ′ of Frost et al. (2022), available for 1995–2020, and the values of F _SG from individual geomagnetic index pairs.

It should be noted that the correction for excess flux from strahl electrons by Frost et al. (2022) provides a slightly larger OSF than the kinematic correction, which was used in past reconstructions by Lockwood et al. (2014b). On the other hand, Figure 11 shows that proper allowance for Parker spiral means we need to reduce OSF by a factor of (1 + μ _ϵ ) ≈ 1.2. These two opposing effects happen to be of similar magnitude, and so, the OSF values derived here are quite similar to those derived by Lockwood et al. (2014b).

The signed OSF derived this way is presented in Figure 12 and compared with the values of Frost et al. (2022) indicated by the blue dots after 1995 in Part A. Before 1995, no suitable strahl data are available, and the blue dots for 1964–1994 are the kinematically corrected values of Lockwood et al. (2009a) that have been adjusted using the best-fit linear regression of the excess flux estimates with the Frost et al. (2022) values for 1995–2020. This composite of satellite OSF values, corrected for excess flux, F S ′ are plotted along the x-axis in Part B. The 2σ uncertainty band of the total OSF reconstruction is slightly broader than that for B and V _SW because it includes the effect of the 5% uncertainty σ _ϵ introduced by the use of the Parker spiral theory on annual timescales (presented in Figure 11D). The correlations for the median of the ensemble for each index pairing are very slightly reduced compared with those for V _SW of B for the same reason. The regression fits of the ensemble medians to the satellite data are similar for the four geomagnetic index pairings in Figure 12B, but those involving IHV give a distinctly greater slope. This is reflected in the reconstruction variations presented in Figure 12A for which the IHV-IDV(1d) and IHV-IDV pairings (the orange and green lines) give the lowest values early in the 20th century. In this context, we note that IHV is not available for cycle 24, and this is a low-activity cycle that helps constrain the reconstructions early in the 20th century when solar activity is low. Even considering this, the agreement is again very good. Supplementary Table S1 in the supplementary material file provides analysis of the correlations of reconstructed B, V _SW , and F _S for each of the four usable index pairings individually and for the median of the ensemble of four million reconstructions.

Here, we have used the aa _H rather than the corrected version of aa employed by Lockwood et al. (2014b), which has introduced some differences. The main effect of this change has been the reduction of discrepancies between the reconstructions for different geomagnetic index pairings as well as the reduction of the overall 2σ uncertainties. However, the use of aa _H has made only minimal difference to the variation of the optimum OSF reconstruction. What has had a much greater effect, and is the main reason why the percentage rise in OSF was somewhat lower than those in the results of Lockwood et al. (2014b), is the evolving improvement to the excess flux calculation. This has resulted in the best estimate of the percentage rise in OSF over the first half of the 20th century changing from 100% in the study by Lockwood et al. (2014b) to the 67% found here. The use of strahl electrons is by far the best method devised (Owens et al., 2017b; Frost et al., 2022), and it generates slightly lower excess flux values than the kinematic correction employed by Lockwood et al. (2014b). This raises all OSF estimates, and a small rise in OSF estimates for early in the 20th century explains the change in the percent rise estimate.

Figures 1 and 2 use running boxcar averages over 11 years to look at the long-term variations. However, because the solar cycles are not exactly 11 years in duration, this gives a residual of approximately 11-year period in the smoothed variations. In addition, the variation in solar cycle durations means that removing them completely using a low-pass filter causes changes to the long-term variation. Lockwood and Frölich (2007) devised a method to give the long-term variations without this ripple nor the need to assume the solar cycle duration. The boxcar running means are taken over a range of interval durations τ between 8 and 14 years. Twice in each solar cycle, these means form “nodes,” where their value is almost the same for all τ. We take the average of the running means at these nodes and interpolate between them using the piecewise cubic hermite interpolating polynomial (PCHIP) interpolation. Figure 13B shows the result for the (revised) international sunspot number, R _ISN . The nodes are shown as yellow points, and the interpolated variation is black. Note that the green line (the boxcar smoothed variation for τ = 11 years) is always very close to the black line, but its residual 11-year ripple has been ironed out. The plot starts shortly after the emergence from the Maunder minimum. Compared with the relatively quiet latest values ( ⟨ R ISN ⟩ τ ≈ 20 ) , we can see three clear maxima: one either side of the Dalton minimum (peaks around 1782 and 1840) and the 20^th-century grand maximum that shows twin peaks (in 1951 and 1985, the former being the larger) and has similar maximum values to the other two but lasts longer. The first of the two peaks in the 20th century grand maximum is higher in ⟨R _ISN ⟩. If we look at all peaks in sunspot number sequences, they are in 1731, 1782, 1840, 1859, 1951, and 1985, giving separations of 51, 58, 19, 92, and 34 years. Hence, it is hard to argue for a regular period beyond the 11-year cycle.

The ISN data are annual mean values. These are available from a full set of daily values for 1818 onward. Before 1818, there are gaps in the daily series, but monthly means are still generated from 1749 onward and annual means from 1700. Hence, the sequence shown is not homogeneous in the data available. In addition, sunspot observers, their apparatus, and their location differ, and indeed, observers’ acuity may vary over time, for example, improving with new instrumentation or decaying with failing eyesight. For later years in the sequence, the mean of observation distribution can be taken, but going back in time, the sequence often relies on single observers, the results from whom need calibration. However, at times, there is little overlap to allow this calibration (Clette and Lefèvre, 2016; Muñoz-Jaramillo and Vaquero, 2019). These calibrations are daisy-chained, and errors in the regression between successive observers (see Lockwood et al., 2016a) accumulate as one goes back in time. A method to calibrate early data against one common standard was devised by Usoskin et al. (2016) using “active day fractions” (ADFs, the fraction of days on which a given observer detects sunspots), which works well at moderate sunspot numbers but has limitations at both low and high levels (Willamo et al., 2018). Chatzistergos et al. (2017) used both daisy-chaining and ADFs to generate a long sunspot number data series, R _GC ; they also used group sunspot numbers, which, unlike R _ISN , means that the more unreliable observations of individual spots are avoided. Figure 13A shows the variation of ⟨R _GC ⟩, and it has similarities—but also some noticeable differences—to R _ISN . In particular, ⟨R _GC ⟩ is increasingly lower than ⟨R _ISN ⟩ before about 1900, which means that for R _GC , the MGSM is larger than the previous maxima in amplitude as well as in integrated solar activity level.

Parts C to H of Figure 13 show the smoothed variation of other indices and reconstructed parameters in the same format as Parts A and B. However, these do no extend back far enough to compare with the peak before the Dalton minima, nor quite to the peak of ⟨R _ISN ⟩ in 1840, but they do show interesting similarities and differences in behavior to both ⟨R _ISN ⟩ and ⟨R _GC ⟩ (shown linearly regressed to the parameter in question in each panel by the light- and dark-gray areas, respectively).

Figure 13C presents the IDV(1d) index. The plot for IDV is essentially identical (not shown), because the agreement of IDV and IDV(1d) is so close (see the bottom panels of Figure 4). The agreement between the long-term variations in ⟨IDV(1d)⟩ (and hence also of ⟨IDV⟩) and the best linear regressions of ⟨R _ISN ⟩ and ⟨R _GC ⟩ is extremely close but is less good for the earlier data. This could potentially indicate that these two geomagnetic sequences are in error at this time or the sunspot data are in error, or both. However, we urge caution before drawing the conclusion that there is an error in one or both of these data sequences. Close inspection shows that the fall in group sunspot numbers after 1985 has been slightly greater than that in sunspot numbers, and hence, it may be a false assumption that the two should agree that closely. Nevertheless, the divergence in the earliest data is considerable. Figure 13C illustrates that the first peak of the 20th-century grand maximum is the larger for IDV(1d) (and so for IDV also) as it is for the sunspot numbers. The earliest IDV(1d) data are more similar to R _ISN than R _GC .

Figure 13D presents the corresponding plot for the sunspot group area ⟨A _G ⟩ composite of data from the Royal Greenwich Observatory (RGO), solar optical observing network, and Kislovodsk, Pulkovo, Debrecen, Kodaikanal, Rome, Catania, and Yunnan observatories, which covers the period between 1874 and 2019 and was compiled and being intercalibrated by Mandal et al. (2020). The group area ⟨A _G ⟩ shows the 20th-century grand maximum and is very similar in amplitude and form to that in the sunspot numbers. The agreement is very good except for the very earliest A _G data, which largely come from RGO.

Figure 13E presents the variation for the aa _H index, which shows great agreement with the aa _HN and aa _HS subindices, with the Potsdam/Seddin/Niemegk k-index composite and the IHV index (see top panels of Figure 4), for which the corresponding plot is therefore almost identical (not shown). However, the aa _H data, even with the Helsinki/St. Petersburg extension, do not quite reach back to the peak in sunspot numbers in 1840, but the 1848 values do imply that this peak is not as pronounced for aa _H as for R _ISN and is more like R _GC . Hence, where IDV(1d) is more like R _ISN at this time, aa _H is more like R _GC , and so using the geomagnetic data to try to discriminate between the two sunspot records gives ambiguous results. However, the most interesting feature is that the second of the 20th-century grand maximum peaks in aa _H is larger than the first, opposite to the sunspot data. Looking at the reconstructed values, the reason is clear. The IMF strength ⟨B⟩ (Figure 13H) and the OSF ⟨F _S ⟩ (Figure 13F) have long-term variations that are very similar to those for the susnpot numbers. In contrast, the solar wind velocity variation ⟨V _SW ⟩ (Figure 13G) is significantly different from those for the sunspot numbers with a broader maximum and dominant peak in 1985. So the dependence of aa _H on V S W 1.76 found in this study explains why the second grand maximum peak is larger than the first for ⟨aa _H ⟩. The smoothed long-term variation ⟨V _SW ⟩ also shows low values at the start of the reconstruction, potentially explaining why the early values indicate that this peak in ⟨aa _H ⟩ is not as pronounced as that in the sunspot number ⟨R _ISN ⟩.

There are several important points to note about these reconstructions. First, the fact that four independent geomagnetic index pairings generate very similar reconstructions of near-Earth IMF strength B, solar wind velocity V _SW , and OSF F _S (demonstrated by Figures 8, 9, 12, respectively) is highly significant. The fact that none of these parameters depends in any way on how and where the geomagnetic data are measured shows that the indices are generally homogeneous in their calibration and accuracy. Note that the indices can differ from each other, but after 1880, at least, the very close similarities of the reconstructions eliminate the possibility that any such differences are an error; in other words, where the geomagnetic indices differ, they should differ, and correcting them to agree with each other would be a retrograde and erroneous step. It is worth noting that to reach this position, we have had to allow for the effect of the geomagnetic latitudes of the stations and how those latitudes have changed because of secular change in the geomagnetic field. This has been done in full for the aa _H and IDV(1d) indices. The IDV index does make allowance for the geomagnetic latitudes of stations but not for the secular change in those geomagnetic latitudes; however, unlike aa _H and IDV(1d) (which use philosophy B in their construction, i.e., daisy chaining of data series that have been corrected for the effects of secular change in the intrinsic field), IDV uses a basket of stations (philosophy A) that make it global in nature. This works when there are sufficient stations in a uniform network across the globe so that secular changes at one station are counterbalanced by changes at others. We note that this becomes increasingly less valid as we go back in time, and the number of available magnetometer stations decreases. Nevertheless, IDV and IDV(1d) agree to a remarkable extent back to approximately 1880. Before then, we have more faith in IDV(1d) because the lack of stations is exposing the weakness of philosophy A when station numbers are low. For the IHV index, Svalgaard and Cliver (2007) showed that the values did not depend on latitude, as long as the stations were well away from the auroral oval (below 55° geomagnetic latitude). Hence, no correction for latitude was needed because stations that were well removed from the auroral oval were chosen at all times.

The agreement between the four reconstructions for B, V _SW , and F _S from the different index pairings shown here means that we can move on from discussions about the accuracy of any of these geomagnetic index data series and exploit the differences between them to infer the behavior and different aspects of the Sun and near-Earth space. The only caveat to this relates to the earliest geomagnetic data. We find very good agreement between the reconstructions after 1880, which leaves a 48-year interval between the advent of the first well-calibrated magnetometer (the work of of Gauss in 1832) and the date when all the indices become unambiguously fit for purpose (i.e., when all the reconstructions based on them agree). For the earliest data, philosophy A (using a basket of stations) is evidently inadequate, and indices using philosophy B (using a well-calibrated and stable station with model allowance for its changing geomagnetic latitude) should be deployed.

Figure 14 presents a comparison of the reconstruction of annual means in OSF presented here with some other recent reconstructions. In Parts A, C, and D, the green and yellow dashed lines indicate the previous reconstruction from geomagnetic data by Lockwood et al. (2014b) (LEA14), the mauve and cyan lines indicate the modeled reconstructions from sunspot numbers by Krivova et al. (2021) (KEA21) and Lockwood and Owens (2014) (L&O14), respectively. The pale-mauve area indicates the reconstruction from the ¹⁴ C cosmogenic isotope by Usoskin et al. (2021) (UEA21) (which extends up to 1899). Figure 14B presents histograms of the annual values from these reconstructions and emphasizes that the MGSM has given unusually high OSF values. The reconstruction of UEA21, based on ¹⁴ C data from 1971 to 1899 gives a mode value considerably lower than that observed for the more recent data. Comparison with the reconstructions based on sunspot number (for after 1616) suggests that there were prior grand minima to 1616 that were deeper than the Maunder minimum. The distribution of reconstructed values presented here (for after 1844) has a higher mode value than for the longer data series derived from ¹⁴ C and sunspot numbers, and for the spacecraft observation period, it is higher still. The distributions for the spacecraft observations, most of which were during the MGSM, contain the highest OSF values in the last 1000 years, and hence, the space age and the MGSM have been an interval of unusually high OSF.

The corrections to the open flux estimates made here (using a better estimate of the excess flux from strahl electron observations and the better allowance for the effects of the Parker spiral gardenhose angle) have important effects. The green and yellow dashed lines in Parts A and C of Figure 14 indicate that the optimum OSF reconstruction presented here has slightly increased OSF values compared with the previous estimates by Lockwood et al. (2014b). There has also been a very slight decrease in the rise since 1900 introduced by using the aa _H index, with its allowance for the secular drift in the geomagnetic latitude of the stations. The LEA14 values always lie within the 2σ uncertainties of the present values. Figure 14A shows that all reconstructions agree closely on the solar cycle and long-term variations of OSF. The percentages of simultaneous values that lie within the estimated uncertainties of the reconstructions presented in this study are given in Part D; they are 83% for the spacecraft observations, 71% for the KEA17 reconstruction from sunspot numbers, 57% for the L&O14 reconstructions from sunspot numbers, and 46% for the UEA21 reconstruction from ¹⁴ C (which ends before the start of the MGSM).

Implementing the improvements detailed in this study, we find the rise in solar activity levels in the 20th-century results in solar cycle means of the OSF F _S increasing from 2.46 × 10¹⁴ Wb in 1906 to 4.10 × 10¹⁴ Wb in 1949, at the peak of the MGSM. Hence, we find that the rise in open flux was by a factor of 67%. It is interesting to note that the peak OSF of the grand maximum derived here is only slightly lower than those estimated by (Lockwood et al., 1999) and by (Lockwood et al., 2014b). However, the value for cycle 14 (around 1906) found here is larger than those estimated in the aforementioned studies, and this lowers the estimated percentage rise in OSF during the MGSM. This is still considerably larger than the 25% rise estimated by Svalgaard and Cliver (2005) and Svalgaard and Cliver (2010), who actually discussed the near-Earth IMF B rather than the OSF F _S but did not acknowledge the differences between these two parameters (Lockwood and Owens, 2011). We here find that the solar cycle means in the IMF B rose from 5.21 in 1901 to 7.63 nT in 1955, indicating a rise of 46%. Lockwood et al. (2006) showed that the low estimate of Svalgaard and Cliver (2005) for the rise in B was because of the adoption of unsafe linear regression procedures.

The variations shown here are otherwise very similar to those derived by Lockwood et al. (2014b), and the variations in the IMF B and solar wind speed V _SW are very similar indeed. Hence, the addition of the solar cycle 24 data, where available, has strengthened the correlations used but not changed these results to any great extent. The slight reduction in the amplitude of the MGSM on OSF was caused by the more rigorous computation of the excess flux using strahl observations.

The variations confirm the similarity of the long-term variations of the IMF and OSF and the sunspot numbers that have been noted before and so give further confidence that OSF modeling based on sunspot numbers (Solanki et al., 2000, 2002; Mackay and Lockwood, 2002; Vieira and Solanki, 2010; Owens and Lockwood, 2012; Lockwood and Owens, 2014; Owens et al., 2017a; Krivova et al., 2021) are not missing a key external factor, although some (small) adjustments of coefficients may be needed to account for the adjusted OSF reconstruction presented here.

Mursula et al. (2017) suggested that the occurrence of low-latitude coronal hole, and hence fast solar wind at Earth and the consequent rise in average solar wind speed, is enhanced during the declining phase of a grand maximum. Figure 14E supports this as average flow speeds are indeed found to be enhanced (relative to the average sunspot numbers) as the open flux declines after the grand maximum. However, the same is true in the rising phase, albeit to a smaller extent. A more detailed study of the enhanced solar wind flow speed will be presented in Paper II (Lockwood et al., 2022).

The reconstructions of the OSF from geomagnetic activity provide a vital dataset for understanding the long-term evolution of the solar atmosphere and the solar magnetic cycle. In particular, like in situ measurements by satellites in near-Earth orbit or in halo orbits around the L1 Lagrange point, the geomagnetic data are measurements of the ecliptic plane and so at low heliographic latitudes (between approximately −7° and +7°) and usually in the streamer belt. However, the results from the Ulysses mission show that they have implications for the solar atmosphere and inner heliosphere at all solar latitudes. This generalization of localized measurements to a global solar value, such as the OSF, may only be valid to first order, and it is probable that the results from Solar Orbiter and Solar Probe will give refinements that can be implemented. In addition, properly joined-up numerical modeling of the solar atmosphere and inner heliosphere offers great promise of improvements.

That being said, tests have shown that the Ulysses result is certainly valid for computing the OSF to first order (Lockwood et al., 2004; Owens et al., 2008; Lockwood and Owens, 2009), and the reconstructions based on geomagnetic data have already been vital. They extend back long enough to allow us to calibrate the very long data series of cosmogenic abundance data in terms of OSF (Usoskin et al., 2021). However, they also raise an interesting dichotomy. The source surface OSF estimates derived from remote photospheric magnetic field observations and numerical modeling are typically a factor of two lower than estimates derived from the in situ satellite and geomagnetic data (Linker et al., 2017; Wallace et al., 2019). Observations of coronal holes also lead to values that are smaller by a similar factor (Lowder et al., 2014, 2017). OSF has been determined using a number of modeling and observational methods. These range from simple PFSS models (Altschuler and Newkirk, 1969; Schatten et al., 1969; Wang and Sheeley, 1992) to complex magnetohydrodynamic (MHD) coronal models (e.g., Lionello et al., 2009). These methods are based on global magnetograms, maps assembled from full-disk observations of the line-of-sight photospheric field. The methods use a number of (different) assumptions and pre-process the input photospheric data in different ways. This can be crucial, for example, Wang and Sheeley (1995) made PFSS and in situ OSF values consistent by adopting a latitude-dependent saturation level for the magnetograph data. However, this is just one of a large number of potential causes of the discrepancy (Stevens et al., 2012). Recent observations by Parker Solar Probe indicate that the discrepancy is not because of heliocentric distance in the heliosphere: modeled and remotely sensed OSF values are also too small than those derived from in situ observations closer to the Sun (at heliocentric distances down to 0.13AU) (Badman et al., 2021).

In Paper II (Lockwood et al., 2022), we look at how these results from geomagnetic observations can be combined with long data series of solar observations to gain a deeper understanding of the long-term change in the solar corona and heliosphere during the rise and fall of the MGSM. The reconstructed data series presented here and the spacecraft interplanetary data are available in the Supplementary Material attached to this study.