Near-Earth heliospheric magnetic ﬁ eld intensity since 1750: 1. Sunspot and geomagnetic reconstructions

We present two separate time series of the near-Earth heliospheric magnetic ﬁ eld strength ( B ) based on geomagnetic data and sunspot number (SSN). The geomagnetic-based B series from 1845 to 2013 is a weighted composite of two series that employ the interdiurnal variability index; this series is highly correlated with in situ spacecraft measurements of B (correlation coef ﬁ cient, r = 0.94; mean square error, MSE = 0.16 nT 2 ). The SSN-based estimate of B , from 1750 to 2013, is a weighted composite of eight time series derived from two separate reconstruction methods applied to four different SSN time series, allowing determination of the uncertainty from both the underlying sunspot records and the B reconstruction methods. The SSN-based composite is highly correlated with direct spacecraft measurements of B and with the composite geomagnetic B time series from 1845 to 2013 ( r = 0.91; MSE=0.24nT 2 ), demonstrating that B can accurately reconstructed by both geomagnetic and sunspot-based methods. The composite sunspot and geomagnetic B time series, with uncertainties, are provided as supporting information.


Introduction
Long-term solar wind variability has potentially important implications for studies of both space [Barnard et al., 2011;Usoskin, 2013] and terrestrial regional and global climate [Gray et al., 2010;Lockwood, 2012;Ineson et al., 2015]. The solar wind has only been observed directly for about 50 years [Neugebauer and Snyder, 1962]. Until recently, high-confidence time series of solar wind and heliospheric magnetic field parameters based on geomagnetic data only extended back to~1900 [Rouillard et al., 2007;Svalgaard and Cliver, 2007b;Lockwood et al., 2009a;Svalgaard and Cliver, 2010]. In this paper we use the sunspot and geomagnetic data sets to consolidate and extend the time series for near-Earth heliospheric magnetic field (HMF) intensity, B, back to 1750. In the companion paper , we compare the two B series derived here (geomagnetic-based and sunspot-based) with the B series based on cosmogenic radionuclide data.
In his landmark paper on the Maunder Minimum, Eddy [1976] conclusively demonstrated that the Sun is a variable star on long time scales. Lockwood et al. [1999] provided further insight to the nature of such long-term change by using geomagnetic activity indices to show that the Sun's open magnetic flux underwent significant-factor of 2-changes during the course of the last century. The Lockwood et al. study reinvigorated the field of long-term solar variability. After~15 years of lively research based on sunspot data [e.g., Solanki et al., 2000Solanki et al., , 2002Wang et al., 2005;Tapping et al., 2007;Balmaceda et al., 2009;Schwadron et al., 2010;Vieira and Solanki, 2010;Vieira et al., 2011;Owens and Lockwood, 2012;Clette et al., 2014; emerging consensus on the reconstruction of B based on geomagnetic data has been forged for the last century [Svalgaard and Cliver, 2010;Lockwood and Owens, 2011]. As noted by Lockwood and Owens [2011], "This is a significant development because individually, each [method] has uncertainties introduced by instrument calibration drifts, limited numbers of observatories, and the strength of the correlations employed." Figure 1, adapted from Lockwood and Owens [2011], shows a comparison of previously published estimates of B from geomagnetic data (described in the figure caption) with the near-Earth spacecraft values. At the time of study, there was reasonable agreement among various geomagnetic reconstructions for the last 100 years, though differences were still present. Subsequent versions have converged yet further. For example, Lockwood et al. [2014a] used four different combinations of geomagnetic indices to better separate the effects of HMF and solar wind speed and place 2σ uncertainties on each: for 1901 (the year with the greatest disagreement), they derived B = 4.7 ± 0.5 nT, in general agreement with the~4.1 nT of Svalgaard and Cliver [2010]. Svalgaard [2014] has replaced Bartel's early diurnal range values with more appropriate data before 1880. Thus, these changes have resolved many of the larger discrepancies that were present (Figure 1), as will be demonstrated in section 2.2.
The present study was initiated as part of an International Space Science Institute (ISSI) international team organized by Leif Svalgaard, Mike Lockwood, and Jürg Beer. It extends the consensus solar wind time series for B back in time to 1750. The present century-long consensus is based primarily on geomagnetic data and, to a lesser extent, on the paleocosmic ray record, discussed in Paper 2. Here we also consider various reconstructions based on the sunspot number for the period from 1750 to the present. The era of consistent highcadence, uninterrupted sunspot coverage began with Schwabe in 1826 while comparable geomagnetic data originated with Gauss, Weber, and the Göttingen Magnetic Union~1830 followed by the British Magnetic program in the early 1840s [Cawood, 1979;Cliver, 1994].
In section 2, we discuss the principal data sets-in situ satellite measurements of B and the sunspot and geomagnetic proxies-and the relationships between the proxies and B. We combine the various time series to obtain the consensus record for each data type and perform the error analysis to determine the uncertainties. The results are summarized and discussed in section 3.

In Situ Spacecraft Observations
The most direct measurement of the near-Earth HMF intensity, B, comes from in situ magnetometer measurements made by spacecraft in the solar wind. The OMNI dataset [King and Papitashvili, 2005] (see also http:// omniweb.gsfc.nasa.gov/) collates the near-Earth solar wind measurements from numerous spacecraft, most recently IMP8, Wind, and ACE. This record begins in 1963 and continues to the present day. Figure 2 (top) shows annual means of B computed from 1 h measurements of the scalar B, referred to as B [OBS] in the remainder of this study. Figure 2 (middle) shows the data coverage for each year. Figure 2 (bottom) shows an estimate of the error in B[OBS] resulting from the incomplete data coverage. This is computed by applying the data gaps for a given year (Y) to all other years in the OMNI interval and recalculating the yearly means.  Lockwood and Owens [2011]. Black dots show in situ spacecraft data from the OMNI near-Earth data set. The yellow line shows the Lockwood et al. [1999] estimate of open solar flux converted to B. The red and blue lines are reconstructions by Lockwood et al. [2009b] and Rouillard et al. [2007]. The grey line and shaded region shows the reconstruction from Svalgaard and Cliver [2010] with uncertainties calculated by Lockwood and Owens [2011].
The error in year Y is then taken to be the average percentage error across all other years. This approach incorporates both the magnitude of the data coverage and the distribution of data gaps. For example, at the time of data access, 2015 had just over 20% data coverage, which is obviously all in the first part of the year (note that while this could be easily updated to use the most recent data, the geomagnetic and sunspot records to which we will compare are currently complete only to the start of 2014). This results in an estimated percentage error higher than for, e.g., 1985, which has similar data coverage but more evenly distributed throughout the year. The effect of the data gaps and averaging timescale on correlations with other parameters was studied by Finch and  who showed that optimum correlation required piecewise removal during data gaps of the data that is to be correlated with B [OBS]. This can be done to a large extent for geomagnetic data because the annual means are computed from higher resolution data. However, the same could not be done for cosmogenic isotope data because the time resolutions achieved are typically 1 year or lower.
Note also that at high time resolutions (seconds to months) geomagnetic activity depends on the southward component of the HMF in the Earth's magnetic frame rather than on the magnitude of the HMF. Stamper et al. [1999] showed that the orientation factor averaged to almost a constant on annual timescales so that the HMF magnitude could be inferred from the geomagnetic data on these timescales. However, even for annual means, this is not quite a constant factor and the small fluctuations place a limit on the maximum correlation between HMF and geomagnetic activity that can be obtained: for annual timescales this limit is approximately 0.95 [Lockwood, 2013].
As estimated errors, particularly after 1965, are relatively small, we compare the various B reconstruction methods to the whole 1963-2013 OMNI data set. For comparison, however, correlations (r) and mean square errors (MSEs) are also computed over the 1995-2013 interval, identified throughout the paper by subscript "95," when spacecraft data coverage is essentially complete.

Geomagnetic Data
Geomagnetic activity provides information about both B and the solar wind flow speed [e.g., Russell, 1975;Svalgaard, 1977;Feynman and Crooker, 1978;Cliver et al., 1998;Lockwood et al., 1999]. Here we use two complementary approaches, both based on the interdiurnal variation (IDV) index [Svalgaard et al., 2003;Cliver, 2005, 2010]. Quoting from Svalgaard and Cliver [2010], "the interdiurnal variability (IDV) index for a given geomagnetic observatory ("station") [is] the average difference without regard to sign, from one day to the next, between hourly mean values of the horizontal component, H, and measured 1 h after midnight [as this is the time most sensitive to solar wind variations]. The average should be taken over a suitably long interval of time, such as 1 year, to eliminate various seasonal complications. IDV has the useful property of being independent of solar wind speed and is highly correlated with [in situ measurements of] the near-Earth heliospheric magnetic field (HMF) strength B." The current version of IDV in this study is from Svalgaard [2014] and is termed S2014. As with earlier versions, it is based on multiple stations, where possible, to increase the signal-to-noise ratio. Figure 6 in Svalgaard [2014] shows that before~1875, there were only four stations or fewer, dropping to a single station (Helsinki) for about a third of the time from 1845 to 1875. After 1950, IDV is based on~20 stations. Combining multiple observatories greatly reduces the "noise" in IDV, particularly during the space age, but this comes at the potential cost of changing characteristics of the time series into the past, as the relative weighting of stations, latitudes, quality, and indices vary with data availability.
The second approach is that recently undertaken by Lockwood and coworkers [Lockwood et al., 2013a[Lockwood et al., , 2013b[Lockwood et al., , 2014a[Lockwood et al., , 2014b, hereafter LEA2013, which uses IDV from daily means, IDV(1d), of H. This definition is equivalent to that for the Bartels' u index [Mayaud, 1980], except that IDV(1d) is based on a single station at any one time rather than on a network of stations and IDV(1d) corrects for the effects of secular change in the geomagnetic latitudes of the stations used. While IDV mitigates the effect of the regular diurnal magnetic variation by considering only data for the midnight hour, IDV(1d) does so by averaging over all 24 h a day. The IDV(1d) index is based on a single sequence of data from geomagnetically similar locations, namely, Helsinki and Eskdalemuir, with limited gaps filled by Potsdam and Seddin. The latter stations are not ideal in that they are at lower latitude than both Helsinki and Eskdalemuir; however, the magnetic latitudes of all these stations (particularly of Eskdalemuir) also change considerably with time because of the secular variation in the geomagnetic field. IDV(1d) makes allowance for this by using the IGRF and gufm1 models of the main field (to estimate the geomagnetic latitude of each station in any one year) along with a study of the effect of geomagnetic latitude on the u values in recent decades (when many stations are available) and normalizes the whole data sequence to the geomagnetic latitude of Eskdalemuir in the year 2000. Thus, IDV(1d), unlike either IDV or u indices, makes allowance for both the secular variation in the geomagnetic field and the effects of changing the distribution of sites used in compiling the index. The advantage of this approach is that the characteristics of the time series are relatively constant with time. The potential disadvantage is that it is prone to inhomogeneities at a single station [e.g., Lockwood et al., 2014b;Svalgaard, 2014] if the data are not properly vetted, though this also applies to IDV for a few intervals prior to~1875 when this sequence also relies on just one or two stations. Because it does not use all of the available data in compiling the index, IDV(1d) can use other data as independent tests. Although no one station can provide a test of the whole period, various intervals are covered by different stations and the agreement is good in all cases (with the exception of Greenwich before 1880 for which the required temperature corrections to the magnetometer data cannot satisfactorily be made).
Comparing and combining the B estimates based on LEA2013 and S2014, denoted B[LEA2013] and B[S2014], has the potential to draw on the individual strengths of the two approaches. As we shall see, the two indices IDV and IDV(Id) are highly correlated. For both IDV and IDV(1d) the correlation coefficient with B[OBS] is very close to the maximum possible (≈0.95) [Lockwood et al., 2014a].  Table 1 also give the correlation coefficients (r) and mean square error (MSE) over the whole 1963-2013 interval, as well as over the 1995-2013 interval, indicated by the subscript "95." The shaded areas show the uncertainty bands, determined as follows: We find the best geometric-mean fit [York et al., 2004]

between B[LEA2013] (or B[S2014]) and B[OBS]
, then introduce a random error to each point, drawn from a two-sided Gaussian distribution, before refitting the data. This is performed 10,000 times and the 90% confidence interval is taken. The best estimate of B is then taken to be the median of the 10,000 fits. This Monte Carlo method of uncertainty determination produces a similar result to the standard geometric-mean fit uncertainty [York et al., 2004], but it more readily allows multiple time series to be combined to produce a weighted ensemble or composite, as discussed below.
The top lines of Table 1 list the statistical properties of the best B[LEA2013] and B[S2014] time series relative to the OMNI data over the 1963-2013 interval. (Note that here we work exclusively with whole-year averages of all data, while [Lockwood et al., 2013a[Lockwood et al., , 2013b computed correlations using piecewise removal of data gaps from the OMNI and geomagnetic data and reported a correlation of 0.95 over the space age. The latter method is preferable for correctly characterizing B because otherwise the geomagnetic annual means include the influence of intervals which have no effect on the HMF data because they were not measured. In this study, however, the use of whole-year averages is necessary in order to treat all time series in the same manner; we have no means of producing piecewise removal of data gaps for the S2014 series and the sunspot-based estimates.) In general, the S2014 approach produces a better match to the whole-year Journal of Geophysical Research: Space Physics 10.1002/2016JA022529  averages of OMNI data than LEA2013 in terms of r and MSE. As discussed earlier, this is because the former averages modern geomagnetic data over many stations, thereby suppressing measurement noise. However, the problem is that the regression for this calibration data may not be as accurate for earlier times when few (in some intervals only one) stations are available. In addition, as described above, data gaps in B[OBS] will have had some effect. The LEA2013 method also produces highly accurate reconstructions over the space age, but is designed to be homogeneous over the whole interval so it reduces concerns about the effect of the changing distribution of stations as we go back in time and is not influenced by data gaps in B [OBS]. The close agreement between the two is very powerful as the S2014 approach gives correlations almost as high as they possibly could have been (considering the HMF orientation factor) and the LEA approach confirms that neither the changing availability of observing stations nor the data gaps in B [OBS] have greatly influenced the results. The green line and green-shaded area in Figure 3c.  This linear scaling has no effect on the calculation of B from the sunspot records. Note that both R I and R G are inherently data composites, comprising different observers, using different equipment in different locations and following different procedures. Figure 4a shows annual R I (blue) and R G (red) back to 1750. These two time series (and, as a result, the time series for B derived from them) differ significantly before~1885. Thus, sunspot number, be it international or group, can become a free parameter in studies of solar and solar terrestrial physics, with the outcome depending on which indicator of solar activity is used.
A series of community meetings has been organized to review and, where appropriate, revise the sunspot number records [Cliver et al., 2013;Clette et al., 2014;. As a result of this process, discontinuities in the international SSN in~1946 [Svalgaard, 2010[Svalgaard, , 2012 [Svalgaard, 2010;Svalgaard, 2013;Clette et al., 2014;Cliver and Ling, 2016], have been identified. In addition, a post-1980 drift has been detected in R I that has been traced to the reference station in Locarno, Switzerland [Clette et al., 2014. Leussu et al. [2013] have reported an inhomogeneity in R I at 1848 when Wolf replaced Schwabe as the primary observer for this time series. In the present study, we use the published "corrected" series, which are summarized here: An R I time series proposed by Lockwood et al. [2014c], here denoted R IL , employs a 12% upward correction of R I before 1946 for the Waldmeier discontinuity and a downward adjustment of 20% before 1849 as suggested by Leussu et al. [2013]. Note that more recent analysis puts the Waldmeier correction at 12% + 3%, though a simple scaling does not account for observed nonlinear effects [Lockwood et al., 2016a]. In the context of the SSN workshops, Clette et al. [2014] reported a preliminary revision of R I (here denoted R IC ). The R IC series used here  contains a stronger adjustment for the Waldmeier discontinuity (~18%) than that of Lockwood et al. [2016a], an alternative correction to that suggested by Leussu et al. [2013] for the discontinuity in 1848 (multiplying the original R I values from 1849 to 1863 by 1.14) and a post-1980 correction related to the transition of the sunspot curatorship from Zürich to Brussels. It is available from the Sunspot Index and Long-term Solar Observations (SILSO) website (www.sidc. be/silso) as S N , V2.0 of the yearly sunspot record. In order to normalize R IC to R I over 1963-1995, R IC is multiplied by 0.72. R IL does not require normalization, as it is identical to R I over this interval. A new derivation of GSN was produced by Svalgaard and Schatten [2016], here denoted R GS . The general approach Svalgaard and Schatten used was to reduce the number of "daisy chains" by identifying a limited number of contiguous principal observers. Secondary observers were scaled against these principal observers to extend the time span covered by each "backbone," thus allowing them to be compared to each other during periods of overlap. The five backbones (designated Staudach, 1739(designated Staudach, -1822Schwabe, 1794Schwabe, -1883Wolfer-Schmidt, 1841-1944Koyama, 1921Koyama, -2000andLocarno, 1950-2015) for the 1750-present interval are pieced together using linear regression, as well as some rather ad hoc (though likely necessary) measures to bridge the data poor interval between the Staudach and Schwabe backbones. Another independent derivation of GSN was recently produced by Usoskin et al. [2016] based on a calibration of solar observers to the reference dataset of RGO data for the period 1900-1976 using the active-day statistics. This series (here designated R GU ) is completely free of daisy-chaining, in contrast to earlier series. These two new group series, R GU and R GS , are scaled by factors 12.15 and 13.88 in order to match the magnitude of R I over the 1963-1995 interval. Figure 4b shows 11 year running means of the six sunspot records in order to highlight the differences in the long-term trends. For R GS and R IC , (white and pink, respectively) the 1750-1800 period exhibits similar or higher sunspot numbers to the space age, whereas R IL and R GU (cyan and black, respectively) show higher values in the space age than 1750-1800. A new construction of the international/Wolf sunspot number [Friedli, 2016] shows long-term trends closer to R GU than R GS and R IC . The differences between the new records will need to be resolved and may necessitate future revisions of SSN-based B series obtained here. At present, we take the approach that by using multiple sunspot records, we can at least constrain the time-dependent uncertainty in SSN-based estimates of B.

Computing B From Sunspot Records
The next step is to compute B from sunspot records. Two different B reconstruction methods are used. The first, denoted B[SSN ½ ], assumes B is linearly correlated to the square root of sunspot number [Svalgaard and Cliver, 2005;Wang et al., 2005]. We note that regressions between SSN ½ and B produce significantly different best fit parameters at different time scales. In particular, the regression of SSN ½ and B for annual means captures the general variation over the solar cycle [Svalgaard and Cliver, 2005], but it does not capture the long-term cycle-to-cycle variation seen for solar cycle means [Wang et al., 2005]. Thus, it is necessary to produce a "double regression" over both important time scales, such that B[SSN ½ ] = A < SSN ½ > SC + B (SSN ½ À < SSN ½ > SC ) + C, where SSN ½ is the square root of the annual mean sunspot number, <SSN ½ > SC is the square root of the solar cycle mean sunspot number, using the solar cycle start/end dates defined by Owens et al. [2011b] on the basis of the average latitude of sunspots. A, B and C are linear best fit coefficients. This purely empirical approach makes no explicit assumptions about the physical mechanism(s) by which B varies. It is single valued, in that a given sunspot number gives a single value of B, regardless of time history. Best linear fits between B[OBS] and SSN ½ over the period 1963-2013 give Note that R G and R GU are not considered here as they are presently only available up to the end of 1995. Clearly, the solar cycle trend will be poorly constrained over these 4.5 cycles. This issue and the implications for a "floor" in B [Svalgaard and Cliver, 2007a] are addressed in section 2.3.4.
The second method derives B from an estimate of open solar flux (OSF), which is modeled as a continuity equation, where the OSF source rate, S, is assumed to vary in phase with SSN [Solanki et al., 2000;Vieira and Solanki, 2010]. This approach incorporates the effect of time history but is based on a number of assumptions about the physical processes by which OSF is produced and lost. We use the recent model of Owens and Lockwood [2012], hereafter OL2012, which uses observed empirical relations and theoretical constraints to reduce the model to a single free parameter. In brief, this model assumes OSF generation associated with coronal mass ejections (CMEs) [Owens and Crooker, 2006] and hence S varies with CME rate. Prior to the advent of coronagraph observations, sunspot number is used as a proxy for CME rate. Using the best fit between observed CME rate [e.g., Yashiro et al., 2004] and R I gives S = ф (0.234 R I 0.540 -0.00153) Wb per Carrington rotation, where ф = 10 12 Wb, the average closed flux carried by a CME [Lynch et al., 2005;Owens, 2008]. From theoretical predictions, the OSF loss term (L) is assumed to vary with the heliospheric current sheet  Figure 5 shows annual time series of B[OBS] (black) and the sunspot-based reconstructions using the SSN ½ (blue) and OL2012 (red) methods. The 90% confidence intervals are shown as shaded areas, computed using the same method as section 2.2. Only R GS and R IL are shown, as R I and R IC are very similar to those of R IL over the space age. Figure 5 (first to fourth panels) clearly show that both B reconstruction methods, using both group-and international-based records, track the observed solar cycle variation and the trends in cycle-to-cycle variability in B, though the full range of the solar cycle variation is generally underestimated. It is not clear if this is an effect of large hemispheric differences in sunspot number, which are particularly apparent around solar maximum [e.g., Lockwood et al., 2012, Figure 7]. Table 1 [Svalgaard and Schatten, 2016] and (third and fourth panels) R IL [Lockwood et al., 2014c]. R GU -  and R IC -  based reconstructions, not shown, are very similar. (fifth panel) A weighted composite time series created from an ensemble of the six B time series resulting from the two B reconstruction methods applied to the three corrected sunspot number series (i.e., R GS , R IL , and R IC ). equivalent. It is important to note, however, that while the OL2012 approach provides a slightly better match to the space age observations than the SSN ½ method, this does not necessarily mean it will continue to provide a better reconstruction outside this period. We note that all sunspot-based reconstructions struggle to reproduce the B[OBS] variation in cycle 20 as accurately as for other cycles. Figure 2 shows that this cycle has more B[OBS] data gaps than cycles 23 and 24 but fewer than cycles 21 and 22. Data gaps or errors in calibration of early B[OBS] data do not appear to be adequate explanations of the cycle 20 differences. Furthermore, geomagnetic reconstructions, which do not suffer from data outages in this period, do match B [OBS] in this cycle (Figure 3).

Journal of Geophysical Research: Space Physics
We now construct a composite B[SSN] time series based on the three corrected SSN series with complete coverage of the space age (R GS , R IL , R IC ). As with the geomagnetic reconstructions, the individual time series are weighted by their "track record," namely, the inverse of the mean square error relative to B[OBS], as listed in Table 1. Figure 5  Similarly, the composite MSE is lower than any individual value. Thus, here the ensemble approach both improves the "best" estimate and reduces the sensitivity to any one method, meaning it is likely to be more robust outside the "training interval" of 1963-2013 than any one individual method (this is further demonstrated in section 2.3.4). More importantly, it provides an estimate of when the different reconstructions do/do not agree, allowing us to estimate uncertainty in the reconstructions resulting from both the methods and the sunspot records. Note, however, that the training interval for the models is the same as the validation interval, which is not desirable. In this instance, this approach is the best possible with the available data, as subdividing the OMNI data into a training and validation intervals leaves insufficient data to adequately constrain the models and Figure 2 shows that the effect of data gaps would be considerably higher in one half of the B[OBS] data series than the other.

Comparing B[SSN] With B[GEO]
While B[SSN], the optimum sunspot-based reconstruction of B, captures over 80% of the observed variation in B over the OMNI era, differences between the various corrected sunspot records before this time mean it can not necessarily be expected to provide the same level of accuracy prior to 1963. This same caveat does not apply to the geomagnetic reconstructions; therefore, we use the B[GEO] composite to recalibrate the weighting factors for the B[SSN] composite; essentially, we are assuming B[GEO] provides a "ground truth" for these purposes. In particular, this is desirable for constraining the cycle-to-cycle variation in the B [SSN ½ ] method.
The best fits of B[SSN ½ ] to the B[GEO] composite are very similar to those obtained by with comparison with OMNI data. Taking the interval 1845-1996, which is covered by all sunspot records, we obtain: Note that by extending the reconstructions back through the Dalton minimum (~1800-1825), they are being used in a parameter regime slightly outside that of the training interval . The lowest possible B value implied by these empirical relations (i.e., the B achieved when both the solar cycle and annual means of SSN are zero) implies a "floor" in B of approximately 2-3 nT. This is well below the B "floor" value of~4.6 nT estimated by Svalgaard and Cliver [2007a] on the basis of the available B reconstructions at the time. The annual mean of the observed B for 2009 was 3.9 nT, demonstrating the importance of considering both the solar cycle and the cycle-to-cycle variations in B, as well as using regressions with caution outside of the observed parameter regime. That said, there is a~2.8 nT estimate of the floor [Cliver and Ling, 2011;Cliver and von Steiger, 2015] that falls within the range of nominal values given above.  Looking first at the statistical differences between the sunspot records listed in

Summary
During the last~60 years, ground-and space-based observations have provided detailed information regarding the spatial and temporal variability of the solar magnetic fields, providing insight into the operation of the solar dynamo [Charbonneau, 2010;Hathaway, 2010]. In addition, satellites have provided observations of the heliospheric magnetic fields from the vicinity of Earth, to the outer edge of the heliopause, both in and out of the ecliptic [Owens and Forsyth, 2013]. In this paper, we have compared the in situ observations of solar wind with geomagnetic observations and sunspot time series to infer the solar wind  [Lockwood et al., 2013b;Svalgaard, 2014]. Both are based on the interdiurnal variation observed at geomagnetic observatories since 1845 but differ in the number of geomagnetic observatories and data utilized.  [Svalgaard and Schatten, 2016;Usoskin et al., 2016] and international sunspot number [Lockwood et al., 2014c;). A weighted composite of the eight B estimates, using the four corrected sunspot records and both B[SSN] estimation techniques, is likely to provide the most robust estimate B and its uncertainty outside of the development period . Compared with B[GEO], this provides a linear correlation coefficient of 0.91 and mean square error of 0.24 nT 2 , an improvement over any one individual sunspot-based B estimate. This composite B[SSN] record over the period 1750-2013, with uncertainties, is provided in the supporting information.

Discussion
A brief review of research on the solar wind in the past is helpful to put the above results in perspective. In their pioneering paper in 1978, 16 years after the discovery of the solar wind, Feynman and Crooker [1978] used a relationship between the aa geomagnetic index and the product of the solar wind speed (V) and the southward component of the solar wind field strength (B s ) to bound V between 237-395 km s À1 and B s between 0.65 and 1.82 nT for the solar minimum year of 1901. More than twenty years would pass before the seminal paper by Lockwood et al. [1999] that separated the V and B dependencies of the aa index and used the discovery by the Ulysses satellite of the constancy of the radial field with latitude to determine the open magnetic flux of the Sun on an annual basis from 1868 to 1996. The pace of research accelerated. Svalgaard et al. [2003] introduced the IDV index which "has the useful property that its yearly averages are highly correlated with the solar wind magnetic field strength (B) and are independent of solar wind speed (V)". Partly because open solar flux estimates were wrongly interpreted as being linearly proportional to B and partly because the 1901 value of B was erroneously low in some reconstructions, a vigorous debate broke out over the reconstruction of B during the past 100 years Svalgaard and Cliver, 2006]. The debate was largely settled by 2010 when a consensus emerged [Svalgaard and Cliver, 2010;Lockwood and Owens, 2011].

10.1002/2016JA022529
While the composite time series for B from 1845 to 2013 based on geomagnetic data and from 1750 to 2013 based on sunspot data presented here represent significant progress, much remains to be done. With the recent publication of the revised sunspot series of Svalgaard and Schatten [2016] and Lockwood et al. [2014c], the SSN-based B series can be extended back to 1610, although the uncertainties will be unavoidably large for the Maunder Minimum (~1650-1710) period Vaquero et al., 2015]. Moreover, as noted in section 2.3.1, differences between the sunspot number time series proposed by Svalgaard and Schatten [2016] and  with those of Lockwood et al. [2014c] and Usoskin et al. [2016] will need to be resolved and may necessitate a revision of the SSN-based B series presented here. Any such changes will not affect our geomagnetic-based B series. An improvement to B[GEO] would be a proper allowance for the effect of data gaps in the B[OBS] data series, as although that was achieved in the B[LEA2013] series it was not exploited in the present paper. An update/extension of the solar wind speed time series [Rouillard et al., 2007;Svalgaard and Cliver, 2007a] based on the new geo-based B series is in preparation. Geomagnetic data is sporadically available after 1722, making possible further potential extensions of geomagnetic-based B and solar wind speed series in the past, though assessing the quality of such data and intercalibrating them with the current records will be challenging. In Part 2 of this study, the series derived here from geomagnetic and sunspot records are compared to cosmogenic radionuclide records .