Medium-term Predictions of F10.7 and F30 cm Solar Radio Flux with the Adaptive Kalman Filter

Forecasting of solar activity with different lead times from hours to decades is very important for many space weather applications (Pesnell 2012, 2016). In particular, predictions of the solar radio flux at wavelengths of F10.7 and F30 cm are required for planning and performing space operations, re-entry predictions (Bastida Virgili et al. 2017a), collision avoidance (Merz et al. 2017), and computation of orbital lifetime. The evolution of the orbital space debris environment is dependent on the decay rates of the objects over time (Bastida Virgili et al. 2014).

F10.7 cm is the flux density of solar radio emissions at a wavelength 10.7 cm averaged over an hour (Tapping 2013), also called the F10.7 solar radio flux or F10.7 index. It comprises a time-varying mix of different emission mechanisms over the solar disk. Together with the F30 index, it serves as a solar proxy of the ultraviolet (UV) solar emission, which heats the Earth's upper atmosphere (Dudok de Wit et al. 2014; Vourlidas & Bruinsma 2018). In comparison to measurements of the true solar UV irradiance, the F10.7 index has long records. The related changes in the thermospheric density cause variations in the atmospheric drag acting on orbiting objects and need to be accurately modeled via atmospheric models. All of them, including the ISO-recommended models, DTM ⁶ (Bruinsma 2015), NRLMSISE ⁷ (Picone et al. 2002), and GOST ⁸ (GOST 2004), considered and assessed by Bastida Virgili et al. (2017b), require F10.7 or F30 as a prime for the solar input. As was shown in the study by Dudok de Wit et al. (2014), in application to DTM2012, ⁹ F30 showed superior performance for modeling of the thermosphere due to a strong presence of thermal bremsstrahlung in the F30 cm flux, which has a strong counterpart in the UV bands.

Forecasting of solar activity is relevant on different timescales, which are dominated by different processes. Short-term forecast is most often used to refer to forecast lead times of several days (up to 30 days; i.e., up to one solar rotation). Medium-term forecast is a forecast for a period when a trend may emerge in the dynamics of the process (such as due to the recurrence of sunspots every 27 days and related activity due to solar rotation), which can be detected and predicted. Medium-term forecast of solar activity usually means predictions one to several months in advance. Long-term forecast is made over a decade or more, i.e., refers to timescales of solar cycles. Such forecasts typically deal with the prediction of the very general characteristics of a solar cycle such as its amplitude, peak time, and period. Medium- and long-term prediction of continuous monthly solar activity indices, which can be quite variable, is quite a delicate task and the predictions can only refer to some trends in the data. Thus, it is also very important first to smooth the data to better detect the trend, which can be predicted.

Various reviews of methods of predicting solar activity have been published over the last decades (e.g., Hathaway et al. 1999; Hathaway 2009, 2010; Pesnell 2012, 2016), and the most recent and extensive one is by Petrovay (2020). Most of these reviews concentrate on techniques for predicting sunspot number. The current approaches can be divided, in general, into three groups. The first group comprises so-called precursor methods based on the identification of some properties of the current cycle, which have predictive power for the next cycle. Long-term predictions are generally made using precursors such as geomagnetic activity (Ohl & Ohl 1979; Feynman 1982; Gonzalez & Schatten 1988; Thompson 1993; Wilson et al. 1998), polar magnetic fields (Schatten et al. 1978, 1996; Schatten & Sofia 1987; Schatten 2005; Svalgaard et al. 2005; Wang & Sheeley 2009; Muñoz-Jaramillo et al. 2012), and sunspot number records or other solar activity indices (Ramaswamy 1977; Lantos 2006; Cameron & Schüssler 2008; Kane 2008; Podladchikova et al. 2008, 2017; Podladchikova & Van der Linden 2011). The second group combines physics-based predictions using models of surface flux transport and dynamo theory (Nandy & Choudhuri 2002; Dikpati & Gilman 2006; Cameron & Schüssler 2007, 2015; Choudhuri et al. 2007; Henney et al. 2012). Finally, the third group represents a variety of statistical methods, which do not require any knowledge of the physics involved (Macpherson et al. 1995; Fessant et al. 1996; Zhang 1996; Conway 1998; Sello 2001; Aguirre et al. 2008; Brajša et al. 2009; Liu et al. 2010; Kakad et al. 2020) and can be applied for medium- and long-term forecasting. Methods of daily short-term prediction of F10.7 and F30 indices use either modifications of (auto-)regression models (e.g., Dmitriev et al. 1999; Liu et al. 2010; Henney et al. 2012) or neural network techniques (e.g., Huang et al. 2009; Yaya et al. 2017). At the same time, the manual forecast method that is currently used at the Solar Influences Data Center (SIDC, Belgium) for short-term F10.7 prediction may be more accurate (Devos et al. 2014).

Despite the large variety of methods of predicting solar activity discussed above, the operational medium-term prediction of solar radio fluxes (i.e., predictions several months in advance) is available only with SOLMAG, the model for predicting solar and geomagnetic activity employed by ESA's Space Debris Office (Mugellesi-Dow et al. 1993; Bastida Virgili et al. 2014). It is based on the well-known regression technique of the medium-term solar cycle prediction, which was introduced by McNish & Lincoln (1949). It uses averaging over past cycles to form a mean cycle as a first approximation, and the difference between the mean cycle and smoothed sunspot numbers to predict future solar activity. Stewart & Ostrow (1970) introduced a modification of the McNish–Lincoln method to predict the monthly mean values of sunspot numbers. Application of the McNish–Lincoln method to the F10.7 index and geomagnetic indices was presented in Holland & Vaughan (1984) and Niehuss et al. (1996). The authors offered to construct the mean cycle by using a Lagrangian linear regression statistical technique. This modification is also used by SOLMAG.

Podladchikova & Van der Linden (2012) presented one of the latest modifications of the McNish–Lincoln method with the development of an adaptive Kalman filter applied to the prediction of sunspot numbers 1–12 months ahead. The developed algorithm was also implemented to improve the medium-term predictions provided by the Standard and Combined methods. The Standard Method is based on an interpolation of Waldmeier's standard curves and provides the predictions of the evolution of the cycle from the average shape of past solar cycles with a specific cycle maximum (Waldmeier 1968). The Combined Method is based on combination of Waldmeier's idea of standard curves and a prediction of the next cycle peak using the aa geomagnetic index (Denkmayr & Cugnon 1997; Hanslmeier et al. 1999). The improved Kalman filter predictions for all three methods are updated monthly at the Sunspot Index and Long-term Solar Observations (SILSO, http://sidc.be/silso), which is the World Data Center for the production, preservation, and dissemination of the international sunspot number.

In this study, we are focusing on a medium-term prediction of F10.7 and F30 indices by developing a new method called RESONANCE (Radio Emissions from the Sun: ONline ANalytical Computer-aided Estimator). First we develop a method of predicting smoothed F10.7 and F30 indices from 1 to 24 months ahead. The method is based on the further development of the adaptive Kalman filter, which is applied to improve medium-term predictions originally produced by the McNish–Lincoln method. Second, we evaluate the performance of the newly developed prediction technique and compare its results with the current state-of-the-art SOLMAG model for medium-term F10.7 predictions. Finally, to demonstrate one of the main applications of predictions of the solar radio index, we re-evaluate past ESA re-entry campaigns for 602 payloads and rocket bodies, and 2344 objects of space debris that re-entered from 2006 to 2019 June, using the forecasts of F10.7 index from the newly developed method and comparing them with those from the SOLMAG model used by ESA.

We develop a method to predict the F10.7 and F30 indices with lead times of 1–24 months. The prediction algorithm has three main steps. First, to process the short-term fluctuations in the measured radio flux we use a 13-month optimized running mean (described in Section 2.2). Second, we use the McNish–Lincoln method to provide initial predictions of the F10.7 and F30 indices (described in Section 2.1). However, the 13-month running mean of the radio flux, which is used as the input to the McNish–Lincoln method, is available at the price of a delay of 6 months with respect to the current time. The radio flux activity over the last 6 months up to the current time is considered to be unknown because of the noise component that is intrinsic to the available monthly mean data. To remove this drawback, in a third step, we develop an adaptive Kalman filter with identification of noise statistics, which assimilates the monthly mean radio flux data over the last 6 months and provides an estimation of the smoothed radio flux at the current time (Section 2.3). The reconstructed radio flux at the current time becomes a new starting point for the predictions, which allows the McNish–Lincoln method to perform without a 6 months delay and thus produce an improved forecast. Throughout the text, when we use the phrase "initial McNish–Lincoln predictions," we refer to the values directly produced by the McNish–Lincoln method; "McNish–Lincoln + Kalman filter predictions" (M&L+KF) is used to mean the improved predicted values produced by applying the combination of the McNish–Lincoln method and the Kalman filter technique (RESONANCE). We also derive analytical expressions for the errors of both the initial McNish–Lincoln predictions and the McNish–Lincoln + Kalman filter predictions. Since those are not yet available in the literature, we include the error derivation in the Appendix.

2.1. McNish–Lincoln Method

The McNish–Lincoln method is based on the construction of a mean cycle by employing all past solar cycles, here as characterized by the F10.7 and F30 indices. Cycles are aligned on the month of the activity minimum, and the average value is calculated for each time step. The predicted solar radio flux F_m for the future month m is calculated according to the following equation:

$$\begin{eqnarray}&&{F}_{m}={\bar{F}}_{m}+{k}_{{mi}}({F}_{i}-{\bar{F}}_{i}).\end{eqnarray} \tag{ 1 }$$

Here, the mean value ${\bar{F}}_{m}$ serves as a first estimation of the solar radio flux at any future month m. This estimation is refined by adding a correction term to it, namely the difference between the last known smoothed value F_i at month i and the value of mean cycle ${\bar{F}}_{i}$ at month i. k_mi is a correction coefficient identified from the least-squares method, which is represented by Equation (A22) in the Appendix. The prediction error of the McNish–Lincoln method is given by Equation (A47), which we derive in analytical form in the Appendix.

In the next subsection we describe the procedure of data preparation and the construction of a mean cycle required by the McNish–Lincoln method.

2.2. Data Preparation

For our analysis, we use monthly mean data of the F10.7 index from the Ottawa and Penticton Radio Observatories, and the F30 index from the Toyokawa and Nobeyama facilities, which provide the longest and continuous records of solar radio flux. The F10.7 and F30 indices are given in solar flux units (sfu) with one 1 sfu = 10⁻²² W m⁻² Hz⁻¹. To filter out short-term fluctuations in the data and to isolate the component associated with the long-term behavior of the solar cycle, a 13-month running mean is widely used (Hathaway et al. 1999). This procedure represents a boxcar averaging of the monthly mean data centered on a given month with equal weights for months −5 to +5 and a half weight for the months −6 at the start and +6 at the end. However, as pointed out by Hurst (1970), Hathaway (2010), Muñoz-Jaramillo et al. (2012), and Podladchikova et al. (2017), the traditional 13-month running mean does a poor job of filtering out high-frequency variations. Thus, we process the noisy record in radio flux data using a 13-month optimized running mean. The advantages of this method over the traditional 13-month running mean are demonstrated and discussed in Podladchikova et al. (2017). The technique is based on finding a balance between fidelity to the data and the smoothness of an approximating curve. The smoothed values are derived from the minimization of the functional

$$\begin{eqnarray}&&J=\beta \displaystyle \sum _{i=1}^{n}{\left({\hat{F}}_{i}-{F}_{i}\right)}^{2}+\displaystyle \sum _{i=1}^{n-2}{\left({F}_{i+2}-2{F}_{i+1}+{F}_{i}\right)}^{2}.\end{eqnarray} \tag{ 2 }$$

Here, ${\hat{F}}_{i}$ is the monthly mean radio flux, F_i is the smoothed flux at month i, and β is smoothing constant, which determines how closely the approximating curve fits the data. For our analysis we use a smoothing constant of β = 0.01, which, as was shown in Podladchikova et al. (2017), provides an effective filtration of noise. The rationale behind this approach is to satisfy two intrinsically conflicting criteria that are used to examine the quality of the approximation of measurement data by a smoothed curve. The data fidelity is evaluated by minimizing the sum of the squared deviations between the fit and the data (first term in Equation (2)), and the smoothness is evaluated by minimizing the sum of squared second derivatives of the fit curve (second term in Equation (2)). The resulting smoothed values minimizing the functional J can be obtained from the solution of a system of n normal equations with n unknowns, where n is the number of available monthly mean data. The technique has been further developed to be applied as the 13-month optimized running mean (see details in Podladchikova et al. 2017).

Figure 1 shows the rms error (RMSE) of the radio flux predictions 1–24 months ahead obtained by the McNish–Lincoln method (with the mean cycle employing cycles 8–23, see below on data preparation for construction of the mean cycle). Panel (a) shows the RMSE for the F10.7 predictions and panel (b) for F30. The red line represents the errors for the 13-month optimized running mean and the blue line those for the traditional 13-month running mean. As can be seen from Figure 1, the application of the 13-month optimized running mean reduces the errors of radio flux predictions compared to the traditional 13-month running mean, especially for longer prediction intervals. The maximal improvements in prediction accuracy with the 13-month optimized running mean reach 13% for the 10 months forecast of F10.7 and 19% for the 14 months forecast of F30 in cycle 24.

Zoom In Zoom Out Reset image size

Data for both F10.7 and F30 indices are available from 1951 November, which constitutes only six 11-year solar cycles 19–24, and does thus not allow construction of the mean cycle employing all past cycles. However, this data set can be augmented by reconstructing the radio flux back to the first cycle using the smoothed sunspot numbers since the correlation between these time series is very high. To reconstruct the solar radio flux data, Bastida Virgili et al. (2014) used a quadratic fit for SOLMAG. We use a third-order linear regression and obtain the following relationship between the smoothed radio flux data F_i and the smoothed sunspot numbers R_i for each month i:

$$\begin{eqnarray}&&{F}_{i}={\beta }_{0}+{\beta }_{1}{R}_{i}+{\beta }_{2}{R}_{i}^{2}+{\beta }_{3}{R}_{i}^{3}.\end{eqnarray} \tag{ 3 }$$

Here, the vectors of regression coefficients ${\beta }^{{\rm{F}}10.7}\,={\left|66.1404,0.4572,0.0018,4.4602\times {10}^{-6}\right|}^{T}$ for F10.7 and ${\beta }^{{\rm{F}}30}={\left|41.3547,0.3669,7.6089\times {10}^{-4},-2.5785\times {10}^{-6}\right|}^{T}$ for F30 are determined from the least-squares method.

The F10.7 and F30 indices reconstructed on the basis of Equation (3) together with the smoothed sunspot numbers are shown in Figure 2. The monthly sunspot numbers are from SILSO, including the corrections applied in Clette et al. (2016). The dark and light green lines in panel (b) show the smoothed F10.7 and F30 indices based on the measured data, while magenta and red lines represent the reconstructed data using the smoothed sunspot numbers given in panel (a). The correlation coefficient between the measured and reconstructed F10.7 cm flux is 0.99 and the standard deviation is 5.43 sfu for cycles 19–24. For the F30 cm flux the correlation coefficient is 0.98 and the standard deviation is 5.71 sfu.

Zoom In Zoom Out Reset image size

As can be seen from Figure 2, the dynamics of the solar radio flux reveals a wide variety for the different solar cycles in terms of amplitude, duration, shape, steepness of the rise, etc. This variability is taken into account by the construction of the mean cycle. However, an important question to answer is whether to use all the reconstructed cycles until the first one or to reject early measurements since, according to McNish & Lincoln (1949), cycles prior to no. 8 are subject to a high level of uncertainty due to inconsistent observations. Thus, the inclusion of a different number of cycles influences the shape of the mean cycle used in the McNish–Lincoln method. As a consequence, rejection of a certain number of cycles may improve the accuracy of the prediction for particular cycles, while worsening it for others. Additional uncertainty is related to the errors in the reconstruction of the solar radio flux from the sunspot numbers using Equation (3). As shown in Figure 2, the differences in the shape of the measured solar radio flux (dark and light green lines) and sunspot numbers (blue line) are observed in particular during the periods of solar maxima for cycles 21–23. For instance, while the second peak of sunspot cycle 22 is smaller than the first one, the amplitudes of both peaks are comparable for the radio flux data. For cycle 23, we observe almost similar peaks in the sunspot data, whereas in the radio flux data the first peak is significantly smaller than the second. In light of this, in Section 3.1, we discuss the performance of the proposed algorithm for different mean cycles employing cycles 8–24 and only the ones where radio measurements are available (cycles 19–24).

2.3. Improvement of McNish–Lincoln Predictions with an Adaptive Kalman Filter

The McNish–Lincoln prediction uses the latest available smoothed radio flux at the price of a 6 months delay with respect to the current time. Figure 3 illustrates an example, when the last 13-month running mean is available for 2010 July (solid red line), and this value with the 6 months delay compared to the current time (2011 January) becomes a starting point for the McNish–Lincoln predictions (dashed black line). The method does not directly use the last available six monthly mean flux data (solid blue line) because of the stochastic component intrinsic to them, though they give significant information about the evolution of the radio flux in the future. We propose to improve the initial McNish–Lincoln predictions by taking into account the dynamics of the monthly mean radio flux over the last six months. To achieve this goal, we develop an adaptive Kalman filter with identification of noise statistics to reconstruct the smoothed radio flux curve up to the current time (solid green line). The estimated radio flux at the current time becomes a new starting point for the improved McNish–Lincoln prediction (dashed blue line), which now performs without a 6 months delay.

Zoom In Zoom Out Reset image size

2.3.1. Model Description

The Kalman filter requires that an investigated dynamical system can be described by a mathematical model in state-space representation, which reflects the changes in the system states. We present the state-space model of the solar radio flux activity in the following way:

$$\begin{eqnarray}&&{F}_{i}={{\rm{\Phi }}}_{i,i-1}{F}_{i-1}+{w}_{i},\qquad (i=1,2,\ldots ,6).\end{eqnarray} \tag{ 4 }$$

Here, F_i is a state vector of the smoothed solar radio flux to be estimated over the last 6 months up to the current time. The transition state matrix Φ_i,i−1 relates a current state vector F_i−1 to a one-step-ahead predicted state vector F_i . The determination of the matrix Φ_i,i−1 usually comes from a physical model of the process under study, which in our case is not specified. However, to provide the best estimation of the unknown solar radio flux F_i the identification of matrix Φ_i,i−1 is crucially important. To get the estimates of the matrix Φ_i,i−1, we suggest using the initial McNish–Lincoln predictions ${\hat{F}}_{i}^{\mathrm{init}}$, which reflect the evolution of the solar radio flux in the future:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{i,i-1}=\displaystyle \frac{{\hat{F}}_{i}^{\mathrm{init}}}{{\hat{F}}_{i-1}^{\mathrm{init}}},\qquad (i=2,3,\ldots ,6).\end{eqnarray} \tag{ 5 }$$

The transition state matrix that relates the last available 13-month running mean F₀ to the one-month-ahead prediction ${\hat{F}}_{1}^{\mathrm{init}}$ is estimated as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{1,0}}=\displaystyle \frac{{\hat{F}}_{1}^{\mathrm{init}}}{{F}_{0}}.\end{eqnarray} \tag{ 6 }$$

However, the transition state matrix Φ_i,i−1 constructed in this way may have unknown errors related to the deviations of the initial McNish–Lincoln predictions ${\hat{F}}_{i}^{\mathrm{init}}$ from the smoothed radio flux activity F_i . The presence of model errors arising from the imperfect description of the dynamics of the process can cause significant estimation errors. Therefore, we add to the state Equation (4) the uncorrelated and unbiased model noise w_i with unknown variance ${\sigma }_{{w}_{i}}^{2}$, which describes the random errors of the initial McNish–Lincoln predictions.

Let ${F}_{i}^{m}\,(i=1,2,\ldots ,6)$ denote the monthly mean radio flux data available over the last 6 months up to the current time. To relate the available measurements F_i ^m to the state vector F_i we need to introduce the following measurement equation:

$$\begin{eqnarray}&&{F}_{i}^{m}={F}_{i}+{\eta }_{i},\qquad (i=1,2,\ldots ,6).\end{eqnarray} \tag{ 7 }$$

Here, η_i is an uncorrelated and unbiased measurement noise with unknown variance ${\sigma }_{{\eta }_{i}}^{2}$, which describes the random deviations of monthly mean measurements F_i ^m from the smoothed flux F_i .

Equations (4) and (7) represent a stochastic state-space model of the solar radio flux activity based on the initial McNish–Lincoln predictions. However, the statistical characteristics of model and measurement noise w_i and η_i introduced into the model are unknown. The application of a Kalman filter provides an optimal estimation with respect to the minimization of the covariance matrix of the estimation error independent of the distribution of noise (Kalman 1960; Sage & Melsa 1971). However, the correct specification of model error statistics is very important since uncertainties in the model and measurements may cause large errors in the prediction and lead to failure of a Kalman filter algorithm. Therefore, to provide the best estimation of the smoothed radio flux F_i , we identify the unknown variances ${\sigma }_{{w}_{i}}^{2}$ and ${\sigma }_{{\eta }_{i}}^{2}$ of the model and measurement noise w_i and η_i using the method for consistent identification of noise statistics presented in Podladchikova & Van der Linden (2012). According to Hathaway et al. (1994), the variance of statistical errors in monthly mean sunspot numbers is proportional to their values. We use the same assumption about the radio flux dynamics and consider the variances ${\sigma }_{{w}_{i}}^{2}$ and ${\sigma }_{{\eta }_{i}}^{2}$ to be proportional to the smoothed radio flux F_i :

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{w}_{i}}^{2} & = & {\alpha }_{w}{F}_{i}\\ {\sigma }_{{\eta }_{i}}^{2} & = & {\alpha }_{\eta }{F}_{i}.\end{array}\end{eqnarray} \tag{ 8 }$$

Here, α_w = 0.2 and α_η = 2.6 are the identified coefficients of proportionality. Note that the application of a Kalman filter provides an optimal estimation with respect to the minimization of the covariance matrix of the estimation error independent of the distribution of noise (Kalman 1960; Sage & Melsa 1971). However, we need to know the variances of the model and measurement noise, which we identified taking into account their increase with the increase in the radio flux activity. The obtained estimates of the variable variances ${\sigma }_{{w}_{i}}^{2}$ and ${\sigma }_{{\eta }_{i}}^{2}$, which take into account the heteroscedasticity of the radio flux data, are further incorporated into the Kalman filter algorithm.

2.3.2. Kalman Filter Algorithm

We estimate the smoothed radio flux described by the state-space model (Equations (4) and (7)) using a Kalman filter, which constitutes a very powerful and efficient data assimilation technique to estimate the state of a process (Kalman 1960). The Kalman filter is a recurrent algorithm that performs in two steps:

1. Prediction is performed to estimate the future state vector of smoothed radio flux F_i one step ahead:

$$\begin{eqnarray}&&{\hat{F}}_{i,i-1}={{\rm{\Phi }}}_{i,i-1}{\hat{F}}_{i-1,i-1},\qquad (i=1,2,\ldots ,6).\end{eqnarray} \tag{ 9 }$$

Here, ${\hat{F}}_{i,i-1}$ represents a predicted estimate of state F_i at step i. We use the first subscript i to indicate the time at which the estimation of the smoothed radio flux is made and the second subscript i − 1 to denote the number of monthly mean radio flux observations ${F}_{1}^{m},{F}_{2}^{m},\ldots ,{F}_{i-1}^{m}$ assimilated by the algorithm to make the prediction. ${\hat{F}}_{i-1,i-1}$ represents a filtered estimate of state F_i−1 at the previous time step i − 1. Here, the second subscript i − 1 indicates that the filtered estimate ${\hat{F}}_{i-1,i-1}$ is obtained on the basis of monthly mean radio flux observations ${F}_{1}^{m},{F}_{2}^{m},\ldots ,{F}_{i}^{m}$.

The prediction accuracy of ${\hat{F}}_{i,i-1}$ in the algorithm is characterized by the variance of prediction error ${\sigma }_{i,i-1}^{2}$ in the following way:

$$\begin{eqnarray}&&{\sigma }_{i,i-1}^{2}={{\rm{\Phi }}}_{i,i-1}^{2}{\sigma }_{i-1,i-1}^{2}+{\sigma }_{{w}_{i}}^{2},\qquad (i=1,2,\ldots ,6).\end{eqnarray} \tag{ 10 }$$

Here, ${\sigma }_{i-1,i-1}^{2}$ is the variance of filtration error at step i − 1. To initiate the Kalman filter procedure, we use the last available value of monthly mean radio flux at the current time as an estimate of the initial state ${\hat{F}}_{\mathrm{0,0}}$ with a variance of the prediction error ${\sigma }_{0,0}^{2}=0$.

2. Filtration is done to incorporate a new measurement F_i ^m in order to obtain an improved estimate of state F_i at step i:

$$\begin{eqnarray}&&{\hat{F}}_{i,i}={\hat{F}}_{i,i-1}+{K}_{i}({F}_{i}^{m}-{\hat{F}}_{i,i-1}),\qquad (i=1,2,\ldots ,6).\end{eqnarray} \tag{ 11 }$$

The estimation accuracy of the filtered value ${\hat{F}}_{i,i}$ is characterized by the variance of the filtration error ${\sigma }_{i-1,i-1}^{2}$:

$$\begin{eqnarray}&&{\sigma }_{i,i}^{2}=(1-{K}_{i}){\sigma }_{i,i-1}^{2}.\end{eqnarray} \tag{ 12 }$$

Here, K_i is the filter gain responsible for the relative weight given to the current measurement F_i ^m and the predicted state ${\hat{F}}_{i,i-1}$ defined as

$$\begin{eqnarray}&&{K}_{i}=\displaystyle \frac{{\sigma }_{i,i-1}^{2}}{{\sigma }_{i,i-1}^{2}+{\sigma }_{{\eta }_{i}}^{2}}.\end{eqnarray} \tag{ 13 }$$

The unknown variances of model and measurement noise ${\sigma }_{{w}_{i}}^{2}$ and ${\sigma }_{{\eta }_{i}}^{2}$ in Equations (10) and (13) are identified from Equation (8), where a filtered estimate ${\hat{F}}_{i-1,i-1}$ at step i − 1 is used instead of the unknown radio flux F_i .

The developed adaptive Kalman filter algorithm assimilates the monthly mean radio flux data over the last 6 months and, as a final output, produces the filtered estimate ${\hat{F}}_{\mathrm{6,6}}$, which represents an estimate of the unknown smoothed radio flux F_i at the current time. The RMSEs of this estimate for cycles 19–24 for both F10.7 and F30 indices are shown in Table 1 (the rows "M&L+KF"). To get the estimate of the smoothed radio flux F_i at the current time with the McNish–Lincoln method, we need to produce a 6 months lead forecast, because the last 13-month running mean, which is used as the input to the McNish–Lincoln method, is available only with a 6 months delay with respect to the current time. The corresponding RMSEs of the McNish–Lincoln predictions are shown in the rows "M&L." The column "ΔRMSE" shows the advantages (in percent) of estimating the smoothed solar radio flux using the proposed Kalman filter compared to the initial 6 months forecast with the McNish–Lincoln method. As shown in Table 1, the developed adaptive Kalman filter reduces the estimation errors of the smoothed solar radio flux at the current time by 23%–46% for F10.7 cm and by 29%–49% for F30 cm monthly mean data.

Table 1. Rms Errors

Cycle No.	Prediction Method	F10.7 cm		F30 cm
		RMSE	ΔRMSE	RMSE	ΔRMSE
19	M&L	8.11		5.69
	M&L+KF	5.14	45%	3.14	45%
20	M&L	7.93		4.95
	M&L+KF	4.25	46%	2.53	49%
21	M&L	6.99		3.65
	M&L+KF	4.86	30%	2.6	29%
22	M&L	13.49		7.64
	M&L+KF	7.56	44%	4.22	45%
23	M&L	9.18		6.39
	M&L+KF	5.03	45%	3.32	48%
24	M&L	6.82		5.1
	M&L+KF	5.22	23%	3.1	39%

Note. RMSE (in units of sfu) of the smoothed estimation of the solar radio flux at the current time, which is considered to be unknown, because the last 13-month running mean is available only with a 6 months delay with respect to the current month. The row "M&L" shows the RMSE of the McNish–Lincoln method, for which the estimation of the smoothed radio flux at the current time is obtained by producing a 6 months lead forecast, because the last 13-month running mean, which is used as the input to the McNish–Lincoln method, is available only with a 6 months delay with respect to the current time. The row "M&L+KF" gives the RMSE of the Kalman filter algorithm, which provides the estimation of the smoothed solar radio flux by assimilating the monthly mean radio flux data over the last 6 months. The column "ΔRMSE" shows the improvements (in percent) of estimating the smoothed solar radio flux using the Kalman filter compared to the initial 6 months forecast with the McNish–Lincoln method.

Download table as:  ASCIITypeset image

3.1. Forecast Performance of M&L+KF Compared to the Initial McNish–Lincoln Predictions

The estimation of the smoothed radio flux at the current time using the developed adaptive Kalman filter allows us to effectively remove an essential drawback of initial McNish–Lincoln predictions, when the last 13-month running mean is available only with a 6 months delay with respect to the current time. To produce an improved radio flux forecast 1–24 months ahead, we again apply the McNish–Lincoln method but use the estimated smoothed radio flux at the current time as the starting point for the predictions. We refer to it as McNish–Lincoln + Kalman filter predictions. The prediction error of the McNish–Lincoln + Kalman filter is represented by Equation (A51), and its analytical derivation is given in the Appendix.

Figure 4 shows the RMSE of the radio flux prediction 1–24 months ahead with the McNish–Lincoln + Kalman filter for two different mean cycles used in the McNish–Lincoln method. The red lines give the RMSE for the mean cycle estimated from cycles 8–24, while the blue lines indicate the prediction errors for the mean cycle derived only from the cycles for which radio flux measurements are available (19–24). Note that each cycle that we predict with the developed algorithm is excluded from the construction of the mean cycle. We show the initial McNish–Lincoln predictions with dashed lines (M&L) and the improved ones using the adaptive Kalman filter with solid lines (M&L+KF). The two columns on the left present the RMSE for F10.7, while the two on the right show the prediction errors for F30.

Zoom In Zoom Out Reset image size

As can be seen from Figure 4, in most cases, application of the McNish–Lincoln + Kalman filter (M&L+KF, solid lines) gives a higher prediction accuracy than the initial McNish–Lincoln predictions (M&L, dashed lines). In general, for both methods the prediction error increases with increasing lead time of the prediction. However, while for cycles 22, 23, and 24 (M&L, dashed blue lines) the errors of the initial McNish–Lincoln predictions of F10.7 cm 12–24 months ahead increase steeply, the application of the McNish–Lincoln + Kalman filter (M&L+KF, blue solid lines) provides a significantly higher prediction accuracy. The estimation of the radio flux activity over the last 6 months up to the current time with the adaptive Kalman filter allows us to reduce the uncertainties in the dynamics and thus to increase the prediction accuracy. Note that the different shape of the RMSE as a function of lead time in solar cycle 19 (steep, almost linear increase) compared to the other cycles is probably related to the early radio observations that had less stability.

In some cases, for example in cycle 22, the F10.7 predictions 11–19 months ahead (red lines) and F30 predictions 7–19 months ahead (blue lines) show a better performance with the initial McNish–Lincoln predictions (M&L, dashed lines). We point out that the indicated prediction interval exhibits a decrease in the prediction errors with increasing lead time, which is untypical for forecasting algorithms. A particularity of cycle 22 is that it has two maxima of almost the same value (Figure 5, blue line) separated by a significant Gnevyshev gap (Gnevyshev 1967; Norton & Gallagher 2010). Also, it has the shortest rise from minimum to maximum of any recorded cycle (just 34 months), and after the first peak, we observe a sharp decrease in the radio flux activity. This specific behavior causes a large mismatch between the dynamics of cycle 22 and the shape of the mean cycle (Figure 5, red line), which reaches a maximum right at the Gnevyshev gap. When we predict the radio flux near the Gnevyshev gap (green point 3) starting from the moment indicated by the green point 2 (prediction 9 months ahead), the underestimated value of the mean cycle ${\bar{F}}_{i}$ leads to an overestimation of the difference $({F}_{i}-{\bar{F}}_{i})$ in Equation (1), and thus to an overestimation of the prediction. The application of a Kalman filter, which assimilates the monthly mean radio flux values over last 6 months, may additionally increase the prediction error. The monthly mean data may indicate a future increase in radio flux activity, while it actually decreases due to a Gnevyshev gap. However, if we start the prediction from the moment marked by green point 1 (prediction 17 months ahead), the difference $({F}_{i}-{\bar{F}}_{i})$ is less distorted due to a smaller mismatch between the shape of the mean cycle and the actual activity, which results in a more accurate prediction for larger prediction intervals. A similar untypical decrease of prediction errors with increasing prediction time is also observed for the prediction of F10.7 in cycle 24 (Figure 5, red lines). However, this is not the case for the F30 prediction, because the Gnevyshev gap in cycle 24 is less prominent for F30 than for F10.7.

Zoom In Zoom Out Reset image size

As can be seen from Figure 4, the accuracy of the predictions depends on both the choice of the mean cycle and the peculiarities of a predicted cycle. The mean cycle estimated from cycles 19–24 has the advantage of choosing available measurements only, but has less diversity of the cycles considered (e.g., double-peak solar cycles, long-lasting solar minimum, abrupt rises of solar activity, etc.). The mean cycle constructed from cycles 8–24 gives a better representation of the different cycle shapes and possible amplitudes. In general, the mean cycle constructed from cycles 8–24 provides a higher accuracy of F10.7 and F30 prediction for cycles 19, 20, and 23 (red lines), while the mean cycle estimated from the measured cycles (19–24) shows better performance for cycles 21, 22, and 24 (blue lines). On average, the prediction accuracy is 10% higher for F10.7 and 5% higher for F30 when constructing the mean cycle from cycles 8–24 than from only the measured ones (19–24). The improvements in prediction accuracy with application of the McNish–Lincoln + Kalman filter reach 36% for F10.7 and 39% for F30. In addition to the mean cycle estimated from cycles 8–24 and 19–24, we also tested the performance of the prediction algorithm for the mean cycle derived from all the available cycles 1–24, cycles 12–24, and cycles 14–24. As the final solution, we chose the mean cycle constructed from cycles 8–24, because, on average, it provides better performance of the McNish–Lincoln predictions than other options.

To investigate the accuracy of predictions for different phases of a cycle, instead of dividing the cycle into two phases by the point of the cycle maximum, ascending and descending, we followed the procedure: fit each cycle with a sine function with argument [0: π] and then divide the cycle into four phases for $[0:\tfrac{\pi }{4}]$, $[\tfrac{\pi }{4}:\tfrac{\pi }{2}]$, $[\tfrac{\pi }{2}:\tfrac{3\pi }{4}]$, $[\tfrac{3\pi }{4}:\pi ]$ intervals of the argument. Figure 6 shows the rms error of predictions provided by the M&L+KF method with the mean cycle derived from cycles 8–24 for the four indicated phases of cycles 21–24. The left panels give the RMSE for F10.7, while the right panels show the same for F30. The blue and magenta lines give the RMSE for low phases of the cycle corresponding to the intervals $[0:\tfrac{\pi }{4}]$ and $[\tfrac{3\pi }{4}:\pi ]$. The red and yellow lines show the RMSE for the high phases of the cycle, i.e., the intervals $[\tfrac{\pi }{4}:\tfrac{\pi }{2}]$ and $[\tfrac{\pi }{2}:\tfrac{3\pi }{4}]$. As can be seen from Figure 6, in general the RMSE values are significantly lower in phases of low as opposed to high solar activity. Some exceptions exist for cycles 19, 20, and 21, which are related to their shapes and closeness of the actual radio flux activity to the mean cycle.

Zoom In Zoom Out Reset image size

3.2. Forecast Performance of M&L+KF Compared to SOLMAG Predictions

SOLMAG is a model for prediction of solar and geomagnetic activity employed by ESA's Space Debris Office. It has short-, medium-, and long-term prediction algorithms, where the last two are also based on the McNish–Lincoln method augmented with corrections from Holland & Vaughan (1984). Cycles before 19 are reconstructed based on the sunspot numbers using a quadratic fit (Bastida Virgili et al. 2014).

Figure 7 shows the performance of the developed M&L+KF method in comparison with the SOLMAG method for cycle 24. The red line shows the RMSE of F10.7 predictions for M&L+KF, while the blue line gives the error for SOLMAG. Figure 7 indicates the overall superior performance of the M&L+KF method for all the prediction intervals. The developed approach statistically outperforms the SOLMAG method by 15.5%–66.5% for cycle 24 with a significant enhancement for longer prediction intervals.

Zoom In Zoom Out Reset image size

To show a more detailed comparison of the performance that takes into account the dependence of the prediction accuracy on the phase of the cycle and the prediction interval, we present in Figure 8 a so-called heat map (right panel). The left y-axis of the heat map indicates the month of the predicted F10.7 value, while the x-axis shows how early the forecast is made (1–24 months in advance). A cell in the heat map shows the prediction error of F10.7 for every month of the cycle. In the left panel we plot the time evolution of the 13-month smoothed monthly mean F10.7 index (blue line) together with the predictions 12 months ahead by M&L+KF (red line) and SOLMAG (green line). The shaded area gives the 1σ error, which is derived with Equation (A51) for the M&L+KF predictions, and can be used to provide an uncertainty range for real-time predictions. In the heat map, we calculate the prediction error as the absolute difference between the predicted and actual radio fluxes. Colors corresponding to negative values (blue–turquoise) indicate better performance of the SOLMAG method, whereas colors of positive values (red–yellow) indicate better performance of the developed M&L+KF method.

Zoom In Zoom Out Reset image size

As can be seen from the heat map (Figure 8), the two cycle peaks are in general better predicted with the proposed McNish–Lincoln + Kalman filter. The number of positive cells is around three times larger for the first peak and around five times larger for the second peak than the number of negative ones. However, the prediction accuracy of the Gnevyshev gap interval is mostly better for the SOLMAG method with six times more negative cells than positive ones. The number of positive and negative cells for the ascending and descending phases is comparable for both prediction methods. However, in the early ascending phase M&L+KF clearly outperforms SOLMAG predictions, which show difficulties at the very beginning of a solar cycle.

As can be seen from the left panel of Figure 8, M&L+KF (red line) gives overestimated predictions in the Gnevyshev gap for the reasons discussed in Section 3.1. At the same time, SOLMAG (green line) gives underestimated predictions of the decreasing radio flux after the first cycle peak. Incorporating different numbers of past cycles into the calculation of the mean cycle changes its shape and affects the accuracy of the predictions. For M&L+KF we use all the cycles beginning from cycle 8, while the mean cycle in the SOLMAG method employs all the cycles starting from cycle 12.

In this section, we present a systematic evaluation of re-entry forecast and test the performance of the F10.7 cm predictions on past ESA re-entry campaigns for 602 payloads and rocket bodies, and 2344 objects of space debris that re-entered from 2006 to 2019 June over the full solar cycle. An example of the re-entry prediction campaign for AVUM, the upper stage of a VEGA-01 rocket, coordinated by ESA is presented in Bastida Virgili et al. (2017a).

Re-entry is the return of an object from outer space into and through the gases of the atmosphere. It can be roughly divided into two categories: controlled—when the spacecraft can determine and influence the time and location of re-entry, and uncontrolled—when the decay is driven by natural forces. Most objects are destroyed by atmospheric heating, but surviving fragments of large objects such as heavy science satellites can create a risk to the population and to property on the ground. During 2008–2017 almost 450 large intact objects have re-entered in an uncontrolled way according to Pardini & Anselmo (2019). This corresponds to a mass of about 900 metric tons. Therefore, the evolution of an uncontrolled re-entry over a timespan from several years to a couple of hours of the remaining lifetime has to be monitored and predicted.

However, it is challenging to predict the re-entry epoch because there are many factors and uncertainties affecting the accuracy of predictions. The main difficulties are related to inaccurate tracking of objects with a long period between passes, uncertain evolution of the attitude, calculation of the mass, surface area, and drag coefficient of an object, computation of thermospheric density, and future solar and geomagnetic activity (Pardini & Anselmo 2013). The thermospheric density is calculated via models of the atmosphere (Lemmens et al. 2016; Bastida Virgili et al. 2017b, 2019) that require as an input the extreme ultraviolet solar radiation and parameters of the solar wind, such as velocity, density, and magnetic field. The rate at which these objects re-enter is not constant because atmospheric decay is the main factor that leads to the re-entry of objects with low Earth orbits. The drag acting on objects depends strongly on the density of the thermosphere, which changes in response to the UV irradiance and the geomagnetic activity. At an altitude of 400 km, the density increases by an order of magnitude from solar minimum to solar maximum (Emmert & Picone 2010). These changes cause large variations in atmospheric drag and therefore determine the frequency with which orbiting objects re-enter the atmosphere.

Figure 9 shows the dependence of re-entry frequency on the phase (a) of solar cycle 24 for (b) payloads and rocket bodies and (c) space debris. As can be seen from Figure 9, the number of re-entered objects for both groups is closely related to the level of solar activity—the majority of objects returned in the phase of the cycle maximum. Consequently, from the view of predictions of solar activity, for most of the objects, the accuracy of the re-entry predictions strongly depends on the ability to forecast solar activity in the maximum phase. As shown in Figure 9(c), the re-entry time of space debris closely follows the evolution of the cycle, reacting immediately to the changes in solar activity, whereas payloads and rocket bodies (Figure 9(b)) show also a large number of re-entries in the area of the second maximum and the declining phase of the cycle, which may be related to the time delay of re-entry for large objects.

Zoom In Zoom Out Reset image size

Figure 10 illustrates the procedure of rerunning past ESA re-entry campaigns. Actual measurements of F10.7 index as well as F10.7 predictions by M&K+KF and SOLMAG can be used as an input to the ESA re-entry prediction tool. As an output, we get the predicted epoch of re-entry and compare it to the actual re-entry time.

Zoom In Zoom Out Reset image size

Over the test period 2006–2019, a total of 1635 predictions were made for 602 payloads and rocket bodies, and 9584 predictions for 2344 objects of space debris. The forecast of the re-entry epoch is performed repeatedly for the same object for different lead times. Figure 11 shows how many forecasts are made (as a percentage) depending on the prediction interval for (a) payloads and rocket bodies and (b) space debris. Zero on the x-axis means that the re-entry prediction is made in the same month as the object re-entered. Although the re-entry prediction interval is less than one month, still one needs to produce a 6 months lead prediction of the smoothed radio flux using the McNish–Lincoln method or to reconstruct the smoothed radio flux value at the current time with the Kalman filter. As there are no re-entry predictions with a 24 months lead time, the last bar in panels (a) and (b) corresponds to a 23 months re-entry forecast.

Zoom In Zoom Out Reset image size

Figure 12 shows the distribution of re-entry prediction errors for the case when the re-entry prediction tool uses the actual known F10.7 as input. The top panels (a) and (b) present the error distributions for payloads and rocket bodies, and the bottom panels (c) and (d) those for objects of space debris. The y-axis indicates the number of cases in every category of re-entering objects with a particular prediction error. When the error is positive, the predicted re-entry time is later than the actual re-entry, while negative errors show cases where the predicted re-entry is earlier.

Zoom In Zoom Out Reset image size

The left panels (a) and (c) show the absolute error in days. The right panels (b) and (d) give the relative prediction error, which is normalized by the lead time of the predictions and determined in the following way:

$$\begin{eqnarray}&&{E}_{p}=\displaystyle \frac{{T}_{{\rm{predicted}}{\rm{re}}\text{-}{\rm{entry}}}-{T}_{{\rm{re}}\text{-}{\rm{entry}}}}{{T}_{{\rm{prediction}}}-{T}_{{\rm{re}}\text{-}{\rm{entry}}}}\times 100{\rm{ \% }}.\end{eqnarray} \tag{ 14 }$$

Here, the numerator shows the difference between the epoch of predicted re-entry ${T}_{{\rm{predicted}}{\rm{re}}\text{-}{\rm{entry}}}$ and actual re-entry ${T}_{{\rm{re}}\text{-}{\rm{entry}}}$. The denominator corrects the error by taking into account how far in advance the prediction is made. The larger the interval between the date of the prediction T_prediction and actual re-entry ${T}_{{\rm{re}}\text{-}{\rm{entry}}}$, the less the relative error E_p .

As can be seen from the histograms (Figure 12), the re-entry prediction tool using the actual known F10.7 still provides a large distribution of re-entry prediction errors. Table 2 shows the percentage of re-entry predictions with absolute errors from 0 to 100, 100 to 200, 200 to 365, and beyond 365 days. Additionally, we show the percentage of re-entry predictions with relative errors from 0% to 33%, 33% to 66%, 66% to 100%, and beyond 100%. Most of the re-entry predictions (58%) in the category of payloads and rocket bodies have an error in the range from 0 to 100 days. However, still 16% of cases have an error of re-entry prediction of more than 365 days, and 17% of cases show a prediction error larger than the lead time. In the category of space debris, the majority of re-entry predictions (61%) have an error larger than 365 days, and 70% of cases have prediction error larger than the lead time. Moreover, for this category, we observe a systematic bias in the re-entry predictions: in 87% of cases, the re-entry is observed earlier than its forecast. This might be related to wrong assumptions about masses, cross sections, or other relevant parameters in the drag model that are not well captured for the small debris objects.

Figure 13 shows the re-entry prediction errors for different stages of solar cycle 24. We tentatively divide the cycle into five phases: the ascending phase, a segment around the first peak, the Gnevyshev gap, a segment around the second peak, and the declining phase. In brackets, we show the number of predicted re-entries made for a given phase of the cycle. The left panels show the absolute error (in days) and relative error (as a percentage) for the payloads and rocket bodies. The right panels give the same, but for the space debris. The y-axis gives the percentage of re-entry predictions for each phase of the cycle with a particular prediction error within the ranges indicated in Table 2. Figure 13 confirms the results from Figure 12 and Table 2, and also shows that the re-entry forecast in the category of space debris has an error larger than 365 days or exceeding the lead time mainly during the Gnevyshev gap, the second peak, and the declining phase of the solar cycle. In summary, even perfect knowledge of solar activity does not provide accurate re-entry predictions, and therefore other components of the atmospheric models also need to be improved.

Zoom In Zoom Out Reset image size

To study the influence of the solar activity predictions on the quality of re-entry forecast and to limit the other sources of uncertainty, we use ESA's re-entry tool to model the re-entry epoch with the actual F10.7 index. We further refer to it as "modeled re-entry epoch." We then rerun the re-entry prediction tool but use the predicted F10.7 index as input. Then for each prediction interval, we calculate in how many cases the McNish–Lincoln + Kalman filter provides better re-entry predictions than the SOLMAG method as compared to the modeled re-entry epoch, and vice versa. The fraction of all the forecasts made for each prediction interval showing an advantage of the McNish–Lincoln + Kalman filter (M&L+KF) is presented in Figure 14 in blue, whereas orange shows an advantage of the SOLMAG method. Brown shows as background the fraction of the non-advantageous method, which can be either M&L+KF or SOLMAG, depending on the forecast lead time. The left panels (a) and (c) show the fraction of cases for the modeled re-entry epoch, which allows us to estimate the impact of the F10.7 forecast on the re-entry predictions. For comparison, we also provide the fraction of cases for the actual re-entry time on the right panels (b) and (d), where the predictions are also affected by the unknown deficiencies of the atmospheric models and other input parameters used in the re-entry prediction tool. The top panels (a) and (b) show the results for payloads and rocket bodies, while the bottom panels (c) and (d) are for space debris. Note that the last bar in the left panels (a) and (c) corresponds to a 24 months re-entry forecast, because we have the re-entry predictions also with a 24 months lead time for the modeled re-entry, in contrast to the actual re-entry time.

Zoom In Zoom Out Reset image size

As can be seen from Figure 14(a), the M&L+KF predictions of the F10.7 index used as an input to the re-entry prediction tool provide a larger fraction of more accurately predicted epochs of re-entry in the category of payloads and rocket bodies. The reason for this is that the majority of predicted re-entries (62%) in that category are observed in the segment around the second cycle peak (24%) and the declining phase (38%). This is the area where the M&L+KF in general outperforms the SOLMAG F10.7 predictions (see the heat map in Figure 8). Some cases that show a larger number of more accurate re-entry predictions with the SOLMAG are related to the interval of the Gnevyshev gap (13% of predicted re-entries), and also partially to the ascending phase and end of cycle 23 (19% of predicted re-entries), and the first cycle peak (6% of predicted re-entries).

For the category of space debris, we see an advantage of the M&L+KF mainly for the prediction interval 1–11 months ahead, while SOLMAG is mostly better for longer prediction intervals 12–24 months ahead (Figure 14(c)). The reason for this is that the majority of the predicted re-entries (59%) are observed in the following three parts of the cycle: ascending phase, the area of the first peak, and the Gnevyshev gap between two maxima. While the percentage of all the predicted re-entries for the Gnevyshev gap between two maxima is the same for both categories, it is increased to 34% for the ascending phase and to 12% for the area of the first peak for space debris, which is about twice as large as for the payloads and rocket bodies. At the same time, it is decreased to 22% for the segment around the second peak and to 19% for the declining phase. As shown in the heat map (Figure 8), SOLMAG provides more accurate F10.7 predictions for larger prediction intervals in the late ascending phase and partially around the first cycle peak. This leads to an increase in the number of cases when the re-entry prediction tool, which uses the SOLMAG F10.7 predictions, provides a more accurate forecast of the re-entry epoch. However, as the M&L+KF shows better performance for these phases of the cycle in the case of a shorter prediction interval (1–11 months ahead), the number of more accurate re-entry predictions with the M&L+KF for this interval is higher than for the SOLMAG. Note that the impact of both M&L+KF and SOLMAG radio flux predictions on the accuracy of re-entry forecast depends on the shape of the cycle, and the presence/absence of a significant Gnevyshev gap. The differences in the re-entry predictions for the actual re-entry (right panels (b) and (d) in Figure 14) is related to other components and uncertainties in the atmospheric models, which are not limited to only the solar activity predictions. Superposition of uncertainties in the atmospheric models and the solar radio flux predictions may lead to improvements of re-entry forecast with the SOLMAG method (Figure 14(d)). Further improvement of the re-entry forecast should include both refinements of the atmospheric models and solar activity predictions.

Figure 15 shows examples of re-entry prediction for (a) the second stage of a Delta 3 rocket, (b) Centaur, the second stage of an Atlas 1 rocket, and (c) the third stage of an N-I rocket. The blue line indicates the relative error of the re-entry predictions (Equation (14)), which use the SOLMAG F10.7 predictions as input to the re-entry prediction tool. The red line shows the same for the M&L+KF F10.7 predictions as input.

Zoom In Zoom Out Reset image size

As can be seen from Figure 15(a), the accuracy of re-entry prediction made 23 months before the re-entry of the second stage of the Delta 3 rocket is comparable for both SOLMAG (blue line) and M&L+KF (red line). However, closer to the re-entry, the prediction error related to the M&L+KF strongly decreases (almost to zero), while it is still big for SOLMAG. The modeled re-entry epoch is 2008 December 29. The 5 months lead forecast predicts that the object will re-enter on 2008 December 10 for SOLMAG and on 2008 December 18 for M&L+KF, which leads to an increase in re-entry prediction accuracy by 18 days. The relative prediction error with M&L+KF dropped almost to zero. The error of re-entry prediction for the Centaur (Figure 15(b)) increases with a decrease in the advance time of prediction from 20 to 14 months, and SOLMAG shows better performance than M&L+KF at these intervals. However, closer to re-entry the error of, for instance, the 9 months lead forecast for M&L+KF significantly decreases, while that for SOLMAG still increases. The modeled re-entry epoch for this object is 2014 October 20. According to the 9 months lead forecast, the object will re-enter on 2014 November 22 for SOLMAG, which gives a prediction that is delayed by more than 1 month (the relative prediction error is 11.5%). At the same time, M&L+KF predicts the re-entry to happen on 2014 October 14, which provides an increase in re-entry prediction accuracy by 27 days. The relative prediction error with M&L+KF decreases to −2.1%, which is about five times smaller than that for SOLMAG. The dynamics of the re-entry prediction errors for the third stage of the N-I rocket (Figure 15(c)) is similar to that for the second stage of the Delta 3 rocket (Figure 15(a)). For instance, the 8 months lead forecast shows that the object will re-enter on 2009 February 5 for SOLMAG and on 2009 March 8 for M&L+KF. The modeled re-entry epoch is 2009 March 9. Thus, the use of the M&L+KF provides a re-entry prediction 31 days more accurate than SOLMAG. The relative prediction error with M&L+KF decreases to −0.4%, which is 32 times smaller than that for SOLMAG. These cases reflect the general trend that for longer lead times SOLMAG provides better predictions, whereas for shorter lead times, i.e., nearer re-entry, M&L+KF shows better prediction results (see Figure 14).

The findings and outcomes of this study can be divided into two main groups. First, we developed a prediction technique (M&L+KF) for the F10.7 and F30 indices with lead times of 1–24 months ahead and compared the performance of our approach with the current state-of-the-art SOLMAG prediction method used by ESA. Second, we applied the developed techniques to re-entry predictions by rerunning past ESA re-entry campaigns for payloads, rocket bodies, and space debris that re-entered over the years 2006–2016, covering in total 602 payloads and rocket bodies, and 2344 objects of space debris.

The developed M&L+KF prediction technique of the solar radio flux consists of three main steps. First, to process the short-term fluctuating component, which should be removed to detect the trend, we applied a 13-month optimized running mean technique (Podladchikova et al. 2017). This results in a better forecasting accuracy than the traditional 13-month running mean with a reduction in the prediction errors by 13% for F10.7 and 19% for F30. Second, we used the McNish–Lincoln method to provide initial predictions of the radio flux indices. However, the 13-month running mean of the radio flux, which is used as input to the McNish–Lincoln method, has a 6 months delay with respect to the current time. The radio flux activity over the last 6 months up to the current time is considered to be unknown because of the noise component intrinsic to the available monthly mean data. To remove this drawback, we developed an adaptive Kalman filter with identification of noise statistics, which assimilates the monthly mean radio flux data over the last 6 months and provides an estimation of the smoothed radio flux at the current time. The reconstructed radio flux at the current time becomes a new starting point for the predictions, which allows the McNish–Lincoln (M&L) method to perform without a 6 months delay. We showed that combining the McNish–Lincoln method with the Kalman filter (M&L+KF) improves the prediction accuracy by 36% for F10.7 and 39% for F30 compared to the initial M&L predictions.

It is important to mention that the proposed Kalman filter technique assimilating the monthly mean radio flux data over the last 6 months represents a universal approach, because its application is independent of the original prediction model. This approach can be applied not only to the McNish–Lincoln predictions but also to the radio flux predictions provided by any other forecasting method. It might be of practical interest in the future to develop and test the proposed approach for the Standard and Combined methods, initially designed to predict sunspot numbers. A comparative study of the performance of these methods developed for prediction of solar radio flux would be of practical importance to test the forecasting capacities for cycles with different shapes and behaviors. Potentially, the proposed approach can be applied to any other proxy of solar activity; however, different observables require specific solutions and tuning.

In the current study, the developed M&L+KF method gives an RMSE of predictions in the range 5–27 sfu for F10.7 and 3–16 sfu for F30. This result statistically outperforms the current state-of-the-art SOLMAG prediction method employed by ESA by 15.5%–66.5% and can be recommended for implementation at ESA.

Additionally, we performed a systematic evaluation of the re-entry predictions by rerunning past ESA re-entry campaigns for payloads, rocket bodies, and objects of space debris that re-entered over the last solar cycle. We showed that even using the known F10.7 as input to the re-entry prediction tool, 16% of the cases for payloads and rocket bodies, and 61% of the cases for space debris, have an absolute error of re-entry prediction of more than 365 days. The percentage of cases where the relative prediction error is larger than the lead time is 17% for payloads and rocket bodies and 70% for space debris. Moreover, in the category of space debris, we observe a systematic bias in the re-entry predictions: in 87% of the cases re-entry is observed earlier than its forecast. These results can help to improve the handling of atmospheric model parameters (assumptions of object masses, cross sections, etc.) in the future and to improve the re-entry prediction tool.

We also tested the influence of the solar radio flux predictions on the accuracy of the re-entry forecast, when the predicted F10.7 was used as input to the re-entry prediction tool. We showed that F10.7 predicted by the proposed M&K+KF method provides a larger fraction of more accurate re-entry forecasts in the category of payloads and rocket bodies than the SOLMAG method. For the category of space debris, M&L+KF shows an advantage in the re-entry forecast for prediction intervals 1–11 months ahead, while SOLMAG is mostly better for longer prediction intervals 12–24 months ahead. This is partially related to the accuracy of the F10.7 predictions by M&L+KF and SOLMAG in cycle 24, which have a prominent Gnevyshev gap.

Further improvements of re-entry forecasts need to include both refinements of the atmospheric models (assumptions of object masses, cross sections, etc.) and solar activity predictions. From our study, it is clear that even with "perfect input" on the solar activity (i.e., the measured instead of the predicted F10.7 index), there are still substantial uncertainties in the re-entry predictions, which hints at the other components of the drag model. For improvements in the medium-term solar activity/radio flux predictions, it is very important to have a better forecast of not only the peak but also the shape of a solar cycle. In that respect, a better understanding and prediction of the Gnevyshev gaps is needed, because they determine the activity amplitude/phase during the period of high solar activity, where the drag affecting space objects is strongest.

As was noted by Dudok de Wit et al. (2014), the F30 index can be a better proxy for modeling the thermosphere–ionosphere system than the currently used 10.7 index. At the same time, as was shown in this study, the predication accuracy of F30 itself is higher than that of F10.7 (Figures 4 and 6). This is related to the lower variability of F30 compared to F10.7, with the variance of monthly mean F30 data being almost three times smaller than that of F10.7. Thus, the incorporation of the F30 index into the re-entry prediction tool is expected to improve the forecast of the re-entry epoch.

The authors acknowledge the Nobeyama Radio Observatory (Nobeyama, Japan) and the Dominion Radio Astrophysical observatory (Penticton, Canada) for providing F10.7 and F30 indices, the team of ESA's Space Debris Office for the SOLMAG predictions and data of the past ESA re-entry campaigns, and WDC-SILSO at the Royal Observatory of Belgium (ROB) for the data sets of sunspot numbers. We thank the referee for valuable comments on this study. A.M.V. acknowledges the support of the Austrian Science Fund (FWF): P27292-N20 and I4555.

To derive the errors of the McNish–Lincoln predictions we make the following estimations.

A.1. The Variance of Smoothed Radio Flux Data

Let F_m ^c denote the values of the smoothed radio flux data at any month m of a cycle c. As the mean cycle in the McNish–Lincoln method is constructed starting from cycle 8, c changes from 8 to N_c , where N_c is the number of the last available cycle. We assume that the sequence ${F}_{m}^{c}\,(c=8,9,\ldots ,{N}_{c})$ is uncorrelated. Then the mean value, ${\bar{F}}_{m}$, of smoothed radio flux data at month m over the considered cycles c is estimated as

$$\begin{eqnarray}&&{\bar{F}}_{m}=\displaystyle \frac{1}{{N}_{c}}\displaystyle \sum _{c=8}^{{N}_{c}}{F}_{m}^{c}.\end{eqnarray} \tag{ A1 }$$

The variance of F_m ^c is then represented by

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{m}^{2} & = & \displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}E\left[{\left({F}_{m}^{c}-{\bar{F}}_{m}\right)}^{2}\right]\\ & = & \displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}\left(E\left[{\left({F}_{m}^{c}\right)}^{2}\right]-E\left[2{\bar{F}}_{m}{F}_{m}^{c}\right]+E\left[{\bar{F}}_{m}^{2}\right]\right).\end{array}\end{eqnarray} \tag{ A2 }$$

Here, the symbol E denotes the mathematical expectation. Let us consider every term of Equation (A2) separately. Term 1:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}E\left[{\left({F}_{m}^{c}\right)}^{2}\right] & = & E\left[\displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}{\left({F}_{m}^{c}\right)}^{2}\right]\\ & = & \displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}{\left({F}_{m}^{c}\right)}^{2}.\end{array}\end{eqnarray} \tag{ A3 }$$

Similarly, term 2 can be estimated as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}E\left[2{\bar{F}}_{m}{F}_{m}^{c}\right]=\displaystyle \frac{1}{{N}_{c}-1}2{\bar{F}}_{m}\displaystyle \sum _{c=8}^{{N}_{c}}{F}_{m}^{c}.\end{eqnarray} \tag{ A4 }$$

Taking into account Equation (A1),

$$\begin{eqnarray}&&\displaystyle \sum _{c=8}^{{N}_{c}}{F}_{m}^{c}={N}_{c}{\bar{S}}_{m}.\end{eqnarray} \tag{ A5 }$$

Then, term 2 can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}E\left[2{\bar{F}}_{m}{F}_{m}^{c}\right]=\displaystyle \frac{1}{{N}_{c}-1}2{N}_{c}{\bar{F}}_{m}^{2}.\end{eqnarray} \tag{ A6 }$$

Term 3:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}E\left[{\bar{F}}_{m}^{2}\right]=\displaystyle \frac{1}{{N}_{c}-1}{\bar{F}}_{m}^{2}\displaystyle \sum _{c=8}^{{N}_{c}}1=\displaystyle \frac{1}{{N}_{c}-1}{N}_{c}{\bar{F}}_{m}^{2}.\end{eqnarray} \tag{ A7 }$$

Thus, Equation (A2) representing the variance ${\sigma }_{m}^{2}$ can be rewritten as

$$\begin{eqnarray}&&{\sigma }_{m}^{2}=\displaystyle \frac{1}{{N}_{c}-1}\left[\displaystyle \sum _{c=8}^{{N}_{c}}{\left({F}_{m}^{c}\right)}^{2}-2{N}_{c}{\bar{F}}_{m}^{2}+{N}_{c}{\bar{F}}_{m}^{2}\right]\end{eqnarray} \tag{ A8 }$$

or

$$\begin{eqnarray}&&{\sigma }_{m}^{2}=\displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({F}_{m}^{c}\right)}^{2}-{N}_{c}{\bar{F}}_{m}^{2}}{{N}_{c}-1}.\end{eqnarray} \tag{ A9 }$$

Similarly to Equation (A9), the variance ${\sigma }_{i}^{2}$ is given by

$$\begin{eqnarray}&&{\sigma }_{i}^{2}=\displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({F}_{i}^{c}\right)}^{2}-{N}_{c}{\bar{F}}_{i}^{2}}{{N}_{c}-1}.\end{eqnarray} \tag{ A10 }$$

Here, F_i ^c denote the values of smoothed radio flux data at any month i of a cycle c (c = 8, 9, ..., N_c ), and ${\bar{F}}_{i}$ represents the mean value of smoothed radio flux data at month i over the considered cycles c.

A.2. Construction of Linear Regression

Let us introduce the following residuals denoting the difference between the smoothed radio flux value and the value of the mean cycle at months m and i:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{m}^{c}={F}_{m}^{c}-{\bar{F}}_{m},\qquad (c=8,9,\ldots ,\,{N}_{c})\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}^{c}={F}_{i}^{c}-{\bar{F}}_{i},\qquad (c=8,9,\ldots ,\,{N}_{c}).\end{eqnarray} \tag{ A12 }$$

Then we form the following linear regression equation between ${{\rm{\Delta }}}_{m}^{c}$ and ${{\rm{\Delta }}}_{i}^{c}$:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{m}^{c}={a}_{{mi}}+{k}_{{mi}}{{\rm{\Delta }}}_{i}^{c}+{\varepsilon }_{{mi}}^{c}.\end{eqnarray} \tag{ A13 }$$

Here, ${\varepsilon }_{{mi}}^{c}$ is the uncorrelated noise representing the model error. Let us introduce the vector

$$\begin{eqnarray}&&Y={\left|{{\rm{\Delta }}}_{m}^{8}{{\rm{\Delta }}}_{m}^{9}\cdots {{\rm{\Delta }}}_{m}^{{N}_{c}}\right|}^{T},\end{eqnarray} \tag{ A14 }$$

the vector of estimated parameters

$$\begin{eqnarray}&&{X}_{{mi}}=\left|\begin{array}{c}{a}_{{mi}}\\ {k}_{{mi}}\end{array}\right|,\end{eqnarray} \tag{ A15 }$$

the matrix

$$\begin{eqnarray}H={\left|\begin{array}{ccc}1 & \ldots & 1\\ {{\rm{\Delta }}}_{i}^{8} & \ldots & {{\rm{\Delta }}}_{i}^{{N}_{c}}\end{array}\right|}^{T},\end{eqnarray} \tag{ A16 }$$

and vector

$$\begin{eqnarray}&&\varepsilon ={\left|{\varepsilon }_{{mi}}^{8}{\varepsilon }_{{mi}}^{9}\cdots {\varepsilon }_{m}^{{N}_{c}}\right|}^{T}.\end{eqnarray} \tag{ A17 }$$

Then the regression Equation (A13) can be rewritten as

$$\begin{eqnarray}&&Y={{HX}}_{{mi}}+\varepsilon .\end{eqnarray} \tag{ A18 }$$

Here, ${\sigma }_{\varepsilon }^{2}$ is the variance of noise ε representing the scatter around the regression line.

We determine the estimate of vector X_mi with unknown parameters a_mi and k_mi on the basis of the least-squares method:

$$\begin{eqnarray}&&{\hat{X}}_{{mi}}=\left|\begin{array}{c}{\hat{a}}_{{mi}}\\ {\hat{k}}_{{mi}}\end{array}\right|={\left({H}^{T}H\right)}^{-1}{H}^{T}Y.\end{eqnarray} \tag{ A19 }$$

Then Equation (A19) can be rewritten as

$$\begin{eqnarray}&&{\hat{X}}_{{mi}}=\left|\begin{array}{c}{\hat{a}}_{{mi}}\\ {\hat{k}}_{{mi}}\end{array}\right|=\left|\begin{array}{c}0\\ \displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}\end{array}\right|.\end{eqnarray} \tag{ A20 }$$

Thus, the unknown parameters of regression can be determined as

$$\begin{eqnarray}&&{\hat{a}}_{{mi}}=0\end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray}&&{\hat{k}}_{{mi}}=\displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}.\end{eqnarray} \tag{ A22 }$$

A.3. Variances of Estimation Errors of Coefficients ${\hat{a}}_{{mi}}$ and ${\hat{k}}_{{mi}}$

Let us introduce the matrix

$$\begin{eqnarray}D={\left({H}^{T}H\right)}^{-1}{H}^{T}=\left|\begin{array}{ccc}\displaystyle \frac{1}{{N}_{c}} & \ldots & \displaystyle \frac{1}{{N}_{c}}\\ \displaystyle \frac{{{\rm{\Delta }}}_{i}^{8}}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}} & \ldots & \displaystyle \frac{{{\rm{\Delta }}}_{i}^{{N}_{c}}}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}\end{array}\right|.\end{eqnarray} \tag{ A23 }$$

Then we can write

$$\begin{eqnarray}&&{\hat{X}}_{{mi}}=\left|\begin{array}{c}{\hat{a}}_{{mi}}\\ {\hat{k}}_{{mi}}\end{array}\right|={DY}.\end{eqnarray} \tag{ A24 }$$

The covariance matrix of the estimation error of ${\hat{X}}_{{mi}}$ is given by

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{cov}\left({\hat{X}}_{{mi}}\right)\\ =\,D\mathrm{var}\left(Y\right){D}^{T}={\sigma }_{\varepsilon }^{2}{{DD}}^{T}={\sigma }_{\varepsilon }^{2}\left|\begin{array}{cc}\displaystyle \frac{1}{{N}_{c}} & 0\\ 0 & \displaystyle \frac{1}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}\end{array}\right|.\end{array}\end{eqnarray} \tag{ A25 }$$

The variance of ${{\rm{\Delta }}}_{i}^{c}$ is given by

$$\begin{eqnarray}&&{\sigma }_{i}^{2}=\displaystyle \frac{1}{{N}_{c}-1}\displaystyle \sum _{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}\end{eqnarray} \tag{ A26 }$$

and then

$$\begin{eqnarray}&&\displaystyle \sum _{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}={\sigma }_{i}^{2}\left({N}_{c}-1\right).\end{eqnarray} \tag{ A27 }$$

Thus, taking into account Equation (A27), Equation (A25) can be rewritten as

$$\begin{eqnarray}\mathrm{cov}\left({\hat{X}}_{{mi}}\right)=\left|\begin{array}{cc}\displaystyle \frac{1}{{N}_{c}} & 0\\ 0 & \displaystyle \frac{1}{{\sigma }_{i}^{2}\left({N}_{c}-1\right)}\end{array}\right|.\end{eqnarray} \tag{ A28 }$$

Therefore, the variances of estimates ${\hat{a}}_{{mi}}$ and ${\hat{k}}_{{mi}}$ are presented as

$$\begin{eqnarray}&&\mathrm{var}\left({\hat{a}}_{{mi}}\right)={\sigma }_{\varepsilon }^{2}\displaystyle \frac{1}{{N}_{c}}\end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray}&&\mathrm{var}\left({\hat{k}}_{{mi}}\right)={\sigma }_{\varepsilon }^{2}\displaystyle \frac{1}{{\sigma }_{i}^{2}\left({N}_{c}-1\right)}.\end{eqnarray} \tag{ A30 }$$

A.4. The Variance of Prediction Errors

Let us rewrite Equation (A13) in the following way:

$$\begin{eqnarray}&&{F}_{m}^{c}={\bar{F}}_{m}+{a}_{{mi}}+{k}_{{mi}}\left({F}_{i}^{c}-{\bar{F}}_{i}\right)+{\varepsilon }_{{mi}}^{c}.\end{eqnarray} \tag{ A31 }$$

The forecast to month m from cycle c is then made as

$$\begin{eqnarray}&&{\hat{F}}_{m}^{c}={\bar{F}}_{m}+{\hat{a}}_{{mi}}+{\hat{k}}_{{mi}}\left({F}_{i}^{c}-\bar{{Fi}}\right).\end{eqnarray} \tag{ A32 }$$

The prediction error ${\delta }_{{im}}$ at month m is then determined by subtracting Equation (A32) from Equation (A31)

$$\begin{eqnarray}\begin{array}{rcl}{\delta }_{{im}} & = & {F}_{m}^{c}-{\hat{F}}_{m}^{c}=\left({a}_{{mi}}-{\hat{a}}_{{mi}}\right)+{k}_{{mi}}\left({F}_{i}^{c}-\bar{{F}_{i}}\right)\\ & & +\ {\varepsilon }_{{mi}}-{\hat{k}}_{{mi}}\left({F}_{i}^{c}-\bar{{F}_{i}}\right)\\ & = & {\varepsilon }_{{mi}}+\left({a}_{{mi}}-{\hat{a}}_{{mi}}\right)+\left({k}_{{mi}}-{\hat{k}}_{{mi}}\right)\\ & & \times \,\left({F}_{i}^{c}-\bar{{F}_{i}}\right).\end{array}\end{eqnarray} \tag{ A33 }$$

The variance ${\hat{\sigma }}_{p}^{2}$ of prediction error ${\delta }_{{im}}$ is given by

$$\begin{eqnarray}&&{\hat{\sigma }}_{p}^{2}=\mathrm{var}\left({\hat{\delta }}_{{mi}}\right)={\hat{\sigma }}_{\varepsilon }^{2}+\mathrm{var}\left({\hat{a}}_{{mi}}\right)+\mathrm{var}\left({\hat{k}}_{{mi}}\right){\left({F}_{i}^{c}-\bar{{F}_{i}}\right)}^{2}\end{eqnarray} \tag{ A34 }$$

or

$$\begin{eqnarray}&&{\hat{\sigma }}_{p}^{2}={\hat{\sigma }}_{\varepsilon }^{2}\left(1+\displaystyle \frac{1}{{N}_{c}}+\displaystyle \frac{{\left({F}_{i}^{c}-\bar{{F}_{i}}\right)}^{2}}{{\sigma }_{i}^{2}\left({N}_{c}-1\right)}\right).\end{eqnarray} \tag{ A35 }$$

A.5. The Estimation of Standard Deviation ${\hat{\sigma }}_{\varepsilon }$

As follows from Equation (A13),

$$\begin{eqnarray}&&{\varepsilon }_{{mi}}^{c}={{\rm{\Delta }}}_{m}^{c}-{k}_{{mi}}{{\rm{\Delta }}}_{i}^{c}-{a}_{{mi}}.\end{eqnarray} \tag{ A36 }$$

The estimation of the variance of noise ${\varepsilon }_{{mi}}^{c}$ is presented as

$$\begin{eqnarray}&&{\hat{\sigma }}_{\varepsilon }^{2}=\displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{m}^{c}-{k}_{{mi}}{{\rm{\Delta }}}_{i}^{c}\right)}^{2}}{{N}_{c}-2}.\end{eqnarray} \tag{ A37 }$$

Here we use N_c − 2 in the denominator because we estimate two regression coefficients, a_mi and k_mi . As follows from Equation (A37),

$$\begin{eqnarray}\begin{array}{rcl}{\hat{\sigma }}_{\varepsilon }^{2} & = & \displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}\left({\left({{\rm{\Delta }}}_{m}^{c}\right)}^{2}-2{k}_{{mi}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}+{\left({k}_{{mi}}{{\rm{\Delta }}}_{i}^{c}\right)}^{2}\right)}{{N}_{c}-2}\\ & = & \displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{m}^{c}\right)}^{2}-2{k}_{{mi}}{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}+{k}_{{mi}}^{2}{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}{{N}_{c}-2}.\end{array}\end{eqnarray} \tag{ A38 }$$

Let us consider all the terms in the numerator of Equation (A38). According to Equation (A27), the first term can be presented as

$$\begin{eqnarray}&&\displaystyle \sum _{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{m}^{c}\right)}^{2}={\sigma }_{m}^{2}\left({N}_{c}-1\right).\end{eqnarray} \tag{ A39 }$$

Let us multiply and divide the second term on ${\sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}$. Then

$$\begin{eqnarray}&&2{k}_{{mi}}\displaystyle \sum _{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}=2{k}_{{mi}}\displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}.\end{eqnarray} \tag{ A40 }$$

Taking into account Equations (A22) and (A27), we can rewrite Equation (A40) in the following way:

$$\begin{eqnarray}&&2{k}_{{mi}}\displaystyle \sum _{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}=2{k}_{{mi}}^{2}{\sigma }_{i}^{2}\left({N}_{c}-1\right).\end{eqnarray} \tag{ A41 }$$

As follows from Equation (A41), the variance ${\sigma }_{i}^{2}$ is given by

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{i}^{2} & = & \displaystyle \frac{2{k}_{{mi}}^{2}\left({N}_{c}-1\right)}{2{k}_{{mi}}{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}}=\displaystyle \frac{{k}_{{mi}}\left({N}_{c}-1\right)}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}}\\ & = & \displaystyle \frac{\tfrac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}\left({N}_{c}-1\right)}{{\displaystyle \sum }_{c=8}^{{N}_{c}}{{\rm{\Delta }}}_{m}^{c}{{\rm{\Delta }}}_{i}^{c}}.\end{array}\end{eqnarray} \tag{ A42 }$$

Thus, taking into account Equation (A10),

$$\begin{eqnarray}&&{\sigma }_{i}^{2}=\displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({F}_{i}^{c}\right)}^{2}-{N}_{c}{\left({\bar{F}}_{i}\right)}^{2}}{{N}_{c}-1}=\displaystyle \frac{{\displaystyle \sum }_{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}}{{N}_{c}-1}.\end{eqnarray} \tag{ A43 }$$

Using Equation (A27), we can rewrite the third term of Equation (A38) in the following way:

$$\begin{eqnarray}&&{\left({k}_{{mi}}\right)}^{2}\displaystyle \sum _{c=8}^{{N}_{c}}{\left({{\rm{\Delta }}}_{i}^{c}\right)}^{2}={\left({k}_{{mi}}\right)}^{2}{\sigma }_{i}^{2}\left({N}_{c}-1\right).\end{eqnarray} \tag{ A44 }$$

Thus, Equation (A38) can be rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{\hat{\sigma }}_{\varepsilon }^{2} & = & \displaystyle \frac{\left({\sigma }_{m}^{2}-2{k}_{{mi}}^{2}{\sigma }_{i}^{2}+{k}_{{mi}}^{2}{\sigma }_{i}^{2}\right)\left({N}_{c}-1\right)}{{N}_{c}-2}\\ & = & \displaystyle \frac{\left({\sigma }_{m}^{2}-{k}_{{mi}}^{2}{\sigma }_{i}^{2}\right)\left({N}_{c}-1\right)}{{N}_{c}-2}.\end{array}\end{eqnarray} \tag{ A45 }$$

The standard deviation ${\hat{\sigma }}_{\varepsilon }$ is then determined as

$$\begin{eqnarray}&&{\hat{\sigma }}_{\varepsilon }=\sqrt{\displaystyle \frac{\left({\sigma }_{m}^{2}-{k}_{{mi}}^{2}{\sigma }_{i}^{2}\right)\left({N}_{c}-1\right)}{{N}_{c}-2}}.\end{eqnarray} \tag{ A46 }$$

Taking into account Equation (A35), the standard deviation of the prediction is given by

$$\begin{eqnarray}&&{\hat{\sigma }}_{p}=\sqrt{\displaystyle \frac{\left({\sigma }_{m}^{2}-{k}_{{mi}}^{2}{\sigma }_{i}^{2}\right)\left({N}_{c}-1\right)}{{N}_{c}-2}}\sqrt{1+\displaystyle \frac{1}{{N}_{c}}+\displaystyle \frac{{\left({F}_{i}^{c}-\bar{{F}_{i}}\right)}^{2}}{{\sigma }_{i}^{2}\left({N}_{c}-1\right)}}.\end{eqnarray} \tag{ A47 }$$

Here, k_mi is derived according to Equation (A22).

A.6. Error Derivation for the McNish–Lincoln + Kalman Filter

In this section we derive the error of predictions obtained by the combination of the McNish–Lincoln method with the Kalman filter. To forecast to month m from cycle c we use Equation (A32). However, here we use the estimated ${\hat{F}}_{i}^{c}$ with the Kalman filter instead of the available smoothed radio flux ${F}_{i}^{c}$. We can then present ${F}_{i}^{c}$ as

$$\begin{eqnarray}&&{F}_{i}^{c}={\hat{F}}_{i}^{c}-\zeta .\end{eqnarray} \tag{ A48 }$$

Here ζ is the random noise with variance ${\sigma }_{i,i}^{2}$ derived from Equation (12), which describes the estimation accuracy of smoothed radio flux data ${F}_{i}^{c}$. Rewriting Equation (A33), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\delta }_{{im}} & = & {F}_{m}^{c}-{\hat{F}}_{m}^{c}=\left({a}_{{mi}}-{\hat{a}}_{{mi}}\right)+{k}_{{mi}}\left({\hat{F}}_{i}^{c}-\zeta -\bar{{F}_{i}}\right)\\ & & +\ {\varepsilon }_{{mi}}-{\hat{k}}_{{mi}}\left({\hat{F}}_{i}^{c}-\zeta -\bar{{F}_{i}}\right)\\ & = & {\varepsilon }_{{mi}}+\left({a}_{{mi}}-{\hat{a}}_{{mi}}\right)+\left({k}_{{mi}}-{\hat{k}}_{{mi}}\right)\left({F}_{i}^{c}-\bar{{F}_{i}}\right)-{\hat{k}}_{{mi}}\zeta .\end{array}\end{eqnarray} \tag{ A49 }$$

The variance ${\hat{\sigma }}_{p}^{2}$ of prediction error ${\delta }_{{im}}$ is given by

$$\begin{eqnarray}&&{\hat{\sigma }}_{p}^{2}={\hat{\sigma }}_{\varepsilon }^{2}\left(1+\displaystyle \frac{1}{{N}_{c}}+\displaystyle \frac{{\left({\hat{F}}_{i}^{c}-\bar{{F}_{i}}\right)}^{2}}{{\sigma }_{i}^{2}\left({N}_{c}-1\right)}\right)+{k}_{{mi}}^{2}{\sigma }_{i,i}^{2}.\end{eqnarray} \tag{ A50 }$$

Here, k_mi ² is estimated using Equation (A22) and replaces ${\hat{k}}_{{mi}}^{2}$ for convenience.

Then, taking into account Equations (A46) and (A50), the standard deviation of the McNish–Lincoln + Kalman filter prediction is given by

$$\begin{eqnarray}&&\,{\hat{\sigma }}_{p}=\sqrt{\displaystyle \frac{\left({\sigma }_{m}^{2}-{k}_{{mi}}^{2}{\sigma }_{i}^{2}\right)\left({N}_{c}-1\right)}{{N}_{c}-2}\left(1+\displaystyle \frac{1}{{N}_{c}}+\displaystyle \frac{{\left({\hat{F}}_{i}^{c}-\bar{{F}_{i}}\right)}^{2}}{{\sigma }_{i}^{2}\left({N}_{c}-1\right)}\right)+{k}_{{mi}}^{2}{\sigma }_{i,i}^{2}}.\,\end{eqnarray} \tag{ A51 }$$

Note that Equation (A51) contains an additional term ${k}_{{mi}}^{2}{\sigma }_{i,i}^{2}$ compared to Equation (A47) describing the errors of the initial McNish–Lincoln predictions. However, the initial McNish–Lincoln predictions use the smoothed radio flux delayed 6 months with respect to the current time, so to estimate the accuracy, for instance, of the 1 month lead forecast, we need to calculate ${\hat{\sigma }}_{p}$ using Equation (A47), which corresponds to prediction 7 months ahead. At the same time, the McNish–Lincoln + Kalman filter allows us to make predictions without a 6 months delay, which in general reduces the prediction error ${\hat{\sigma }}_{p}$ calculated with Equation (A51).