Predictability of the Solar Cycle Over One Cycle

The various proxies of solar activity show a roughly 11-year cycle period with widely varying amplitudes. The cycle amplitude gives a rough idea of the frequency of space weather storms of all types, from radio blackouts to geomagnetic storms to radiation storms. For this reason there is a practical need for solar cycle prediction in our technological society. The longer the prediction extends with higher reliability, the more desirable it is for decision makers. Moreover, the predictions, especially the dynamo-based prediction, can also provide us the opportunity to examine our insight into the physical mechanisms underlying the solar cycle.

Existing attempts to predict the level of solar activity can be broadly divided into two groups: (1) those predicting an ongoing cycle, and (2) those predicting the next or future cycle(s). The first group can be further divided into medium-term predictions (months in advance) and long-term predictions (years in advance). Three methods for short-term prediction of ongoing cycles are the McNish–Lincoln method (McNish & Lincoln 1949), the standard method, and the combined method at Sunspot Index and Long-term Solar Observations.³ For long-term ongoing cycle predictions, people usually use the curve-fitting function (Waldmeier 1935; Hathaway et al. 1994; Li 1999; Du 2011; Li et al. 2017) or use the similarity of solar cycles, e.g., Li et al. (2002).

For long-term predictions of the future cycle(s), existing attempts can be broadly divided into empirical methods and methods based on dynamo models. The empirical methods can be further divided into two subgroups. One group is based on extrapolating the past records using a purely mathematical or statistical analysis, which has wide applications in econometrics.⁴ The other subgroup of the empirical methods are the precursor methods, which are based on correlations between certain measured quantities in the declining phase of a cycle and the strength of the next cycle, e.g., Thompson (1993), Svalgaard et al. (2005), and Cameron & Schüssler (2007). Around the end of cycle 23, two groups of people made the dynamo-based solar cycle prediction available for the first time (Choudhuri et al. 2007; Jiang et al. 2007; Dikpati & Gilman 2006; Dikpati et al. 2006). The differently observed data assimilated into prediction models and different understandings of the flux transport mechanisms cause markedly distinct predictions for cycle 24 (Karak et al. 2014). For a review about the solar cycle prediction, see Petrovay (2010) and references therein.

Although there have been numerous attempts to predict the solar cycle, the reliability of the different methods is still controversial. The predictions of cycle 24 using the different prediction methods, with a comparison to the true strength of solar cycle 24, are clearly presented in Pesnell (2012), who analyzed 75 predictions for cycle 24. The proved most-successful method is the precursor method based on the polar field (Schatten et al. 1978; Schatten & Sofia 1987; Layden et al. 1991; Schatten et al. 1996; Svalgaard et al. 2005) and the geomagnetic activity (Ohl & Ohl 1979; Feynman 1982), which is a good proxy for the Sun's polar field during cycle minima (Wang & Sheeley 2009).

The correlation of the polar field at cycle minimum on the subsequent cycle strength is the natural result of the Babcock–Leighton (BL) dynamo (Babcock 1961; Leighton 1969). Recently Cameron & Schüssler (2015) demonstrated that the net toroidal magnetic flux generated by differential rotation within a hemisphere of the convection zone is determined by the emerged magnetic flux at the solar surface. The demonstration and the predictive skill of the polar field indicate that the solar dynamo is of the BL type. In the BL dynamo framework, the poloidal field regeneration, half of the dynamo loop, is accessible to direct observation. Furthermore, the unavoidable time taken to wind up the surface poloidal, means that surface fields observed today will be the source of the sunspot-forming toroidal magnetic component in the future. This justifies the method to use the polar field at cycle minimum to produce a cycle prediction. But can we predict the polar-field evolution before cycle minimum so that the prediction scope can be extended?

The polar-field evolution can be understood in terms of the Surface Flux Transport (SFT) model (Wang et al. 1989; Mackay & Yeates 2012; Jiang et al. 2014b). The tilt angle of the sunspot group, with respect to the E–W direction, leads to each active region contributing to the global axial dipole moment of the Sun and is an essential factor in the polar-field generation. The tilt angle has a systematic component that corresponds to Joy's law and a larger random component, i.e., the tilt-angle scatter. Jiang et al. (2014a) found the observed scatter of the tilt angles causes a variation of 30%–40% in the resulting polar field around activity minima. By including detailed information about the individual tilt angles and magnetic polarities of the bipolar magnetic regions (BMRs) that emerged during cycle 23, the weak polar field around the end of cycle 23 is well reproduced by Jiang et al. (2015). Recently series of BL dynamo models successfully reproduce the variability of the solar cycle by introducing random BMR tilts (Jiang et al. 2013; Olemskoy & Kitchatinov 2013; Cameron & Schüssler 2017; Lemerle & Charbonneau 2017; Nagy et al. 2017; Karak & Miesch 2017). The randomness of the sunspot emergence—part of the toroidal-to-poloidal part of the dynamo loop—makes the dynamo a stochastic process. Even if the poloidal-to-toroidal part of the dynamo loop was fully deterministic, the intrinsic random features of flux emergence limits the scope of the solar cycle prediction. Furthermore, the randomness of the sunspot emergence due to the complex flux emergence dynamics (Weber et al. 2011, 2013) also causes uncertainty, even for the prediction of an ongoing cycle. This randomness raises the question of whether or not long-term predictions with reasonable levels of uncertainty, are possible. A second question is the extent to which we can predict the time dependence of the activity of a solar cycle beyond only predicting the cycle amplitude.

Jiang et al. (2011) show the dependence of the statistical properties of sunspot emergence on the cycle phase and strength based on an empirical analysis of the historical sunspot number data (see also Muñoz-Jaramillo et al. 2015). Random realizations of sunspot group emergences during a cycle can be reconstructed for a given maximum sunspot number based on the statistical properties. The SFT model can be used to produce a large-scale field evolution over the Sun's surface. Cameron et al. (2016) and Jiang & Cao (2018) applied the above scheme in a Monte-Carlo approach to predict the range of the polar field around the end of cycle 24 about 3–4 years before the minimum, including error bars. Similar approaches are also presented in Hathaway & Upton (2016) and Iijima et al. (2017). In this paper, we generalize and update the previous method to investigate the predictability of the solar cycle over one cycle at different phases. Predictions include time evolution of the monthly sunspot number of an ongoing cycle and time evolution of the smoothed sunspot number over one cycle, all with estimated uncertainties. The investigation of the predictability beginning from different starting points distinguishes our scheme from the existing predictions of the solar cycle.

The paper is organized as follows. In Section 2, we study the properties of solar cycle profiles, which will be adopted by Section 3 for the prediction of sunspot emergence for an ongoing cycle. The SFT model, which will be used to describe and to predict the large-scale field evolution over the solar surface, is presented in Section 4. The correlation between the dipole moment at cycle minimum and the subsequent cycle strength based on a homogeneous axial dipole moment data set is given in Section 5. In Section 6, we present the results about the predictability of the subsequent cycle and its application to cycle 25. Our summary and discussion are given in Section 7.

At present there are different time series of the sunspot number, which differ substantially before cycle 12 (Clette et al. 2014). Hence in this study we only investigate the properties of solar cycle profiles during cycles 12–24, and used Sunspot Number Version 2.0.⁵ The data are of the monthly mean total sunspot number (${R}_{\mathrm{mn}}$) and the corresponding 13 months smoothed sunspot numbers (${R}_{\mathrm{sm}}$). The timing of each cycle minimum used throughout the paper is taken from the NGDC.⁶ We first analyze the properties of the smoothed sunspot numbers. The monthly sunspot number will be analyzed at the end of this section.

Figure 1 shows the observed smoothed sunspot number as a function of months from the start of a cycle for cycles 12–24. We define the maximum phase of a cycle as the time period when the sunspot number surpasses 70% of maximum sunspot number of the cycle. The ascending phases, the maximum phases, and the declining phases are denoted in the dotted, dashed, and solid curves in Figure 1, respectively. There are two distinctive features about the shape of the solar cycle. The first is the rising phase, which obeys the Waldmeier effect (Waldmeier 1955), where stronger cycles tend to show a faster rise of activity levels during their ascending phase than weaker cycles (Lantos 2000; Cameron & Schüssler 2008). The second is that once the solar cycle begins to decline, all cycles decline in a similar way (Ivanov & Miletsky 2014; Cameron & Schüssler 2016). A large scatter of the individual cycle about the means over all the cycles is also shown. During the declining phase, the profile of a solar cycle shows noisier short-term variation than the early time period. These properties can be understood under the framework of the BL dynamo. A weaker polar field at the beginning of a cycle generates less of a toroidal field to form the sunspot groups than it would if the polar field was strong. This provides an explanation to the Waldmeier effect and is also demonstrated by Karak & Choudhuri (2011) using a BL dynamo model. Cameron & Schüssler (2016) interpreted the similar decline phases in terms of oppositely directed toroidal flux bands in each hemisphere that diffuse and cancel across the equator when the distance of the center of the activity belts from the equator becomes about equal to their width. Stronger cycles show wider activity belts and thus start to decline earlier than weaker cycles. Nagy et al. (2017) and Kitchatinov et al. (2018) showed that peculiar BMRs with large tilt angles that emerge during the rising phase of a cycle have large effects on the amplitude of the descending phase. Hence the stochastic mechanism due to random emergence of the peculiar BMRs during the rising phase causes the later phases of the cycle to be noisy, as seen in Figure 5.

Zoom In Zoom Out Reset image size

Several of these properties have been exploited in the development of a parameterization of the sunspot number during a cycle in terms of two simple parameters, the starting time (t₀) and the amplitude (a; as given in Hathaway et al. 1994, hereafter HWR94). The determination of the two parameters can be made a few years after the start of a cycle. The function captures the essence of the Waldmeier effect. The rising phase is presented in the form of t³, where t is the time in months from the start of a cycle. We use the function, i.e., Equation (1) of HWR94, to fit the cycles 12–24 individually in the form

$$\begin{eqnarray}&&f(t)=\displaystyle \frac{a{(t-{t}_{0})}^{3}}{\exp [{(t-{t}_{0})}^{2}/{b}^{2}]-c},\end{eqnarray} \tag{ 1 }$$

where $b{(a)=27.12+25.15/(a\times {10}^{3})}^{1/4}$ and c = 0.71. The solar cycle overlap, i.e., new cycle spots appearing at high latitudes while old cycles are still seen at low latitudes, is about 3 years (Wilson et al. 1988; Harvey 1992). We set weights to the first 1.5 years of each cycle by an error function of the form $\tfrac{1}{2}[1+\mathrm{erf}(2\tfrac{t-{t}_{\mathrm{ovp}}}{{t}_{\mathrm{ovp}}})]$, where ${t}_{\mathrm{ovp}}=0.75$ years to decrease the effect of the data during the first 1.5 years on the cycle fit.

Figure 2 shows examples of the cycle fit overplotted with the observed sunspot number evolution from 1878 to the present, based on the fit either 4 years (upper panel) or 8 years (lower panel) after the start of the cycle. Since the function captures the essence of the Waldmeier effect, the fit and smoothed sunspot number are consistent with each other during the ascending phase. We measure the relative error ${s}_{{\rm{f}}2{\rm{o}}}$ between the fit and the observation during the maximum phases for different cycles when the fit is done at different timings t of cycles. Here ${s}_{{\rm{f}}2{\rm{o}}}$ is defined as

$$\begin{eqnarray}&&{s}_{{\rm{f}}2{\rm{o}}}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{\left(\tfrac{{R}_{\mathrm{sm}}({t}_{i})-f({t}_{i})}{f({t}_{i})}\right)}^{2}}{N}},\end{eqnarray} \tag{ 2 }$$

where N is the number of months of the maximum phase for different cycles. The results are shown in Figure 3. The fits to cycles 12–24 are done from 1.5 to 6 years in cycles with 0.5 years interval. The strength of a cycle defined by the maximum sunspot number is indicated by a given color shown by the color bar. The stronger cycles, shown in colors from green to red, tend to have smaller relative errors, which are usually less than 20%, for the fits starting from 2 years onwards. The weaker cycles, shown in from black to cyan, tend to match with the fits more poorly. From 3 years into a cycle, ${s}_{{\rm{f}}2{\rm{o}}}$ is always within 50% relative error. This leads us to only make ongoing cycle prediction from at least 3 years into a cycle throughout the rest of this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

During the later phases, the fitting curves in Figure 2 show larger deviations from the observations. BMRs during the descending phase tend to be located at low latitudes. Jiang et al. (2014a) demonstrated that BMRs with the same initial axial dipole moment at lower latitudes have larger contributions to the polar field. Hence the BMRs during the descending phase play important roles in the polar-field evolution. This motivates us to suggest an improvement to the predictions. We use two methods to measure the differences $\bigtriangleup f(t)$ between the observation ${R}_{\mathrm{sm}}(t)$ and the fitted function f(t). The first is

$$\begin{eqnarray}&&{R}_{\mathrm{sm}}(t)=f(t)+\bigtriangleup f(t),\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&\bigtriangleup f(t)=\overline{\bigtriangleup f(t)}+\bigtriangleup f^{\prime} (t),\end{eqnarray} \tag{ 4 }$$

where $\overline{\bigtriangleup f(t)}$ will be added to function (1) to correct the systematic deviation from the cycle fits especially in the descending phase. Meanwhile, we estimate the fluctuation of the cycle shape by measuring the standard deviation (${\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}$) of fits including the correction using $\overline{\bigtriangleup f(t)}$ from observations. $\overline{\bigtriangleup f(t)}$ and ${\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}$ are calculated by

$$\begin{eqnarray}&&\bar{\bigtriangleup f(t)}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}({R}_{{\rm{sm}}}{(t)}_{i}-f{(t)}_{i})}{N}{\rm{and}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{({R}_{\mathrm{sm}}{(t)}_{i}-f{(t)}_{i}-{\overline{\bigtriangleup f(t)}}_{i})}^{2}}{N}},\end{eqnarray} \tag{ 6 }$$

where the symbol i denotes different cycles 12–24 and N is equal to 13. The random component $\bigtriangleup f^{\prime} (t)$ will be added to predict the sunspot number evolution as a set of random realizations with zero average and the standard deviation of ${\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}$.

The other method to correct the fit is

$$\begin{eqnarray}&&{R}_{\mathrm{sm}}(t)=f(t)+\bigtriangleup r(t)f(t),\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&\bigtriangleup r(t)=\overline{\bigtriangleup r(t)}+\bigtriangleup r^{\prime} (t),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\overline{\bigtriangleup r(t)}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}\left(\tfrac{{R}_{\mathrm{sm}}{(t)}_{i}-f{(t)}_{i}}{f{(t)}_{i}}\right)}{N},\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&{\sigma }_{\bigtriangleup {r}^{{\prime} }(t)}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{\left(\tfrac{{R}_{\mathrm{sm}}{(t)}_{i}-f{(t)}_{i}-{\overline{\bigtriangleup r(t)}}_{i}f{(t)}_{i}}{f{(t)}_{i}+{\overline{\bigtriangleup r(t)}}_{i}f{(t)}_{i}}\right)}^{2}}{N}},\end{eqnarray} \tag{ 10 }$$

where i values are the same as those in Equation (6), which are the 13 cycles from cycle 12 to cycle 24. The random component $\bigtriangleup r^{\prime} (t)$ will be added to predict the sunspot number evolution as a set of random realizations with zero average and the standard deviation of ${\sigma }_{\bigtriangleup {r}^{{\prime} }(t)}$.

Figure 4 shows the time evolution of $\overline{\bigtriangleup f(t)}$, ${\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}$, $\overline{\bigtriangleup r(t)}$, and ${\sigma }_{\bigtriangleup {r}^{{\prime} }(t)}$ when the fits are done from 3 years or more into a cycle. The curves in each panel correspond to the different timings (from the third year to the ninth year into cycles with one year cadence) to do the cycle fits. We see that $\overline{\bigtriangleup r(t)}$ and ${\sigma }_{\bigtriangleup {r}^{{\prime} }(t)}$ increase for the later phase of cycles. When the fits are done at the third year, the fourth year, and the fifth year, $\overline{\bigtriangleup r(t)}$ and ${\sigma }_{\bigtriangleup {r}^{{\prime} }(t)}$ have the almost same profiles. If the fits are done at the sixth year to the eighth year, there is a significant increase of the match between the observed values and the fitting function and the $\bigtriangleup r(t)$-values decrease accordingly. Close to the minimum time period, i.e., from the ninth year into a cycle, there are large deviations among different fits. The small sunspot number during the decline phase contributes to the large variations of the $\bigtriangleup r(t)$-value. In contrast, $\overline{\bigtriangleup f(t)}$ and ${\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}$ during this time period weakly depend on the fits, which are done at different phases of a cycle.

Zoom In Zoom Out Reset image size

The fitting functions for the time evolution of $\overline{\bigtriangleup r(t)}$ and ${\sigma }_{\bigtriangleup r^{\prime} }(t)$ are listed in the following.

$$\begin{eqnarray}&&\overline{\bigtriangleup r(t)}=8.68\exp (-{z}_{1}^{2}/2),t\lt 72\,\mathrm{months},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\overline{\bigtriangleup r(t)}=2.64\exp (-{z}_{2}^{2}/2),72\,\mathrm{months}\leqslant t\leqslant 108\,\mathrm{months},\end{eqnarray} \tag{ 12 }$$

where ${z}_{1}=\tfrac{t-200.44}{45.86}$ and ${z}_{2}=\tfrac{t-149.0}{27.31}$. Whether Equation (11) or Equation (12) is used depends on when the prediction is done. When it is less than 72 months into a cycle, we take Equation (11). When it is more than 72 months and less than 108 months into a cycle, we take Equation (12). The corresponding ${\sigma }_{\bigtriangleup {r}^{{\prime} }(t)}$ has a weak dependence on the fitting time and satisfies the following form:

$$\begin{eqnarray}&&{\sigma }_{\bigtriangleup {r}^{{\prime} }(t)}=30.92\exp (-{z}_{3}^{2}/2),t\leqslant 108\,\mathrm{months},\end{eqnarray} \tag{ 13 }$$

where ${z}_{3}=\tfrac{t-318.73}{77.35}$. During the later phases, $\overline{\bigtriangleup f(t)}$ obeys the following profiles:

$$\begin{eqnarray}\overline{\bigtriangleup f(t)} & = & -2.41\exp (-{z}_{4}^{2}/2)\\ & & \ +\ 36.8-0.2t,t\geqslant 108\,\mathrm{months},\end{eqnarray} \tag{ 14 }$$

where ${z}_{4}=\tfrac{t-119.3}{2.0}$. The corresponding ${\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}$ is

$$\begin{eqnarray}{\sigma }_{\bigtriangleup {f}^{{\prime} }(t)} & = & -3.06\exp (-{z}_{5}^{2}/2)\\ & & \ +\ 47.23-0.314t,t\geqslant 108\,\mathrm{months},\end{eqnarray} \tag{ 15 }$$

where ${z}_{5}=\tfrac{t-108.16}{1.455}$. The two parameters $\overline{\bigtriangleup f(t)}$ and ${\sigma }_{\bigtriangleup {f}^{{\prime} }(t)}$ are shown in the red curves in the lower left and the lower right panels, respectively.

We now analyze the monthly sunspot number ${R}_{\mathrm{mn}}$. The deviation of ${R}_{\mathrm{mn}}$ from the smoothed sunspot number ${R}_{\mathrm{sm}}$ is regarded as a random component. We measure the relative deviations ${s}_{\mathrm{mn}2\mathrm{sm}}$ of ${R}_{\mathrm{mn}}$ from ${R}_{\mathrm{sm}}$. It is defined as

$$\begin{eqnarray}&&{s}_{\mathrm{mn}2\mathrm{sm}}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{\left(\tfrac{{R}_{\mathrm{mn}}{(t)}_{i}-{R}_{\mathrm{sm}}{(t)}_{i}}{{R}_{\mathrm{sm}}{(t)}_{i}}\right)}^{2}}{N}},\end{eqnarray} \tag{ 16 }$$

where i is the 13 cycles from cycle 12 to 24. The result is shown in Figure 5. We see that during the maximum phase, the relative deviations ${s}_{\mathrm{mn}2\mathrm{sm}}$ keep to small values of about 0.2. Then the value increases with time during the decline phase, which means noisier later phases of a solar cycle than earlier phases, the reason for which is given in the second paragraph of this section. During the first 2 years, the large ${s}_{\mathrm{mn}2\mathrm{sm}}$ is due to the cycle overlap. The fitting functions are

$$\begin{eqnarray}&&{s}_{\mathrm{mn}2\mathrm{sm}}=\left\{\begin{array}{l}0.936\exp (-{p}_{1}^{2}/2)+0.196,t\lt 36\,\mathrm{months},\\ 0.375\exp (-{p}_{2}^{2}/2)+0.196,t\geqslant 36\,\mathrm{months},\end{array}\right.\end{eqnarray} \tag{ 17 }$$

where ${p}_{1}=\tfrac{t+12.21}{14.346}$ and ${p}_{2}=\tfrac{t-128.4}{25.357}$.

Zoom In Zoom Out Reset image size

3.1. Prediction of the Sunspot Number Evolution

With the features of the solar cycle shape listed in Section 2, we may predict the time evolution of the monthly sunspot number ${F}_{\mathrm{mn}}(t)$ for the remainder of a cycle when it is some months, denoted as n, entering into the cycle. We separate the sunspot number into the systematic component and the random component. We first use the method described in Section 2 to fit the observed smoothed sunspot number to derive the function f(t). We then get the predicted systematic part of smoothed sunspot number evolution $\overline{{F}_{\mathrm{sm}}}(t)$ = $f(t)+\overline{\bigtriangleup f(t)}$ or $f(t)+\overline{\bigtriangleup r(t)}f(t)$ for the different values of n according to Equations (11), (12), and (14). The random deviations from the systematic component measured by Equations (13) and (15) are determined by using Monte-Carlo simulations. The random realization of the sunspot emergence is denoted as ${F}_{\mathrm{sm}}(t)$, which is equal to $\overline{{F}_{\mathrm{sm}}}(t)+\bigtriangleup f^{\prime} (t)$ or $\overline{{F}_{\mathrm{sm}}}(t)+\bigtriangleup r^{\prime} (t)f(t)$, where $\bigtriangleup f^{\prime} (t)$ and $\bigtriangleup r^{\prime} (t)$ have zero averages and their standard deviations satisfy Equations (13) and (15). The monthly sunspot number ${F}_{\mathrm{mn}}(t)$ is ${F}_{\mathrm{sm}}(t)+{r}_{\mathrm{sm}}(t){F}_{\mathrm{sm}}(t)$, where the standard deviation of ${r}_{\mathrm{sm}}$ satisfies Equation (17) and the average of ${r}_{\mathrm{sm}}$ is zero.

We take two timings, 4 and 8 years, into cycles to compare predictions of ongoing cycles $\overline{{F}_{\mathrm{sm}}}(t)$ and the observations ${R}_{\mathrm{sm}}(t)$. The difference is measured by the goodness-of-fit, which is given by

$$\begin{eqnarray}&&\chi =\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{\left(\tfrac{\overline{{F}_{\mathrm{sm}}}({t}_{i})-{R}_{\mathrm{sm}}({t}_{i})}{\sigma ({t}_{i})}\right)}^{2}}{N}},\end{eqnarray} \tag{ 18 }$$

where $\sigma ({t}_{i})$ corresponds to Equation (6) and N corresponds to the numbers of months from the timing of the prediction to the end of the cycle. When χ is equal to 1.0, it indicates that the fitting function passes within one standard deviation of the observed data points. Table 1 shows χ values of the predicted results of the sunspot number evolution during cycles 12–24 at 4 and 8 years into cycles. The χ values are also calculated for the prediction just based on HWR94 method, i.e., fits using Equation (1). We see that χ values are always much smaller than that of HWR94, especially when the predictions are at the earlier phase, i.e., 4 years into a cycle. This demonstrates the improvement of the predictive skill by the current strategy compared to HWR94.

Table 1.  Measurements of the Accuracy of the Predicted Mean Values of the Smoothed Sunspot Number Evolution During Cycles 12–24 at 4 and 8 Years into Cycles by the Goodness-of-fit χ

	Prediction at 4 years into a Cycle		Prediction at 8 years into a Cycle
	Current Method	HWR94	Current Method	HWR94
cy.12	1.43	1.94	1.10	0.40
cy.13	0.73	1.20	0.66	1.71
cy.14	0.80	1.20	1.32	0.48
cy.15	1.02	1.14	0.91	0.44
cy.16	0.76	1.16	0.88	0.37
cy.17	0.79	1.04	0.37	1.39
cy.18	0.61	1.27	0.13	1.07
cy.19	0.83	1.21	0.78	1.87
cy.20	1.83	2.91	1.53	2.85
cy.21	0.65	1.01	0.70	0.64
cy.22	0.61	0.99	0.51	0.84
cy.23	0.84	1.90	0.46	1.59
cy.24	1.18	1.01	⋯	⋯

Note. The current method given by this paper and the profile given by HWR94 are both listed for comparison.

Download table as:  ASCIITypeset image

Another advantage of the method is to quantify the uncertainty of the predictions, which corresponds to the predictability of the sunspot number evolution. Tables 2 and 3 give the evaluation of the predictive skill of the smoothed sunspot number ${F}_{\mathrm{sm}}(t)$ and the monthly sunspot number ${F}_{\mathrm{mn}}(t)$, respectively. The values show the percentages of the observed sunspot numbers, which are out of 1σ, 2σ, and 3σ values of the predicted variations. For some cycles, e.g., cycle 12 and cycle 20, the percentages are larger than 32% and 4.6% for 1σ and 2σ fluctuations, respectively. But all of the observed sunspot numbers are within a 3σ range of the prediction. The predictive skill for monthly sunspot number is worse than that of the smoothed sunspot number.

Table 2.  Quantitative Measurements of the Prediction of the 13 Months Smoothed Monthly Sunspot Number ${F}_{\mathrm{sm}}(t)$

	Prediction at 4 years into a Cycle			Prediction at 8 years into a Cycle
	1σ	2σ	3σ	1σ	2σ	3σ
cy.12	47.5	11.2	0.0	43.7	0.0	0.0
cy.13	0.0	0.0	0.0	12.2	0.0	0.0
cy.14	6.4	0.0	0.0	54.3	0.0	0.0
cy.15	20.8	5.5	0.0	29.1	0.0	0.0
cy.16	14.9	0.0	0.0	15.4	0.0	0.0
cy.17	18.4	0.0	0.0	0.0	0.0	0.0
cy.18	8.2	0.0	0.0	12.0	0.0	0.0
cy.19	17.7	0.0	0.0	4.0	0.0	0.0
cy.20	70.3	37.3	0.0	67.4	13.9	0.0
cy.21	13.3	0.0	0.0	14.8	0.0	0.0
cy.22	5.9	0.0	0.0	0.0	0.0	0.0
cy.23	20.6	0.0	0.0	0.0	0.0	0.0
cy.24	39.5	10.4	0.0	⋯	⋯	⋯

Note. The values show the percentages of the observed sunspot number, which are out of the 1σ, 2σ, and 3σ values of the predicted variations.

Download table as:  ASCIITypeset image

Table 3.  Quantitative Measurements of the Prediction of the Monthly Total Sunspot Number ${F}_{\mathrm{mn}}(t)$

	Prediction at 4 years into a Cycle			Prediction at 8 years into a Cycle
	1σ	2σ	3σ	1σ	2σ	3σ
cy.12	47.5	17.5	8.7	43.7	0.0	0.0
cy.13	29.8	3.1	0.0	40.8	0.0	0.0
cy.14	32.0	16.0	9.3	17.4	0.0	0.0
cy.15	40.3	9.7	5.5	50.0	4.2	4.2
cy.16	39.2	5.4	1.4	34.6	0.0	0.0
cy.17	39.4	5.3	0.0	28.6	3.6	0.0
cy.18	28.8	6.8	0.0	28.0	0.0	0.0
cy.19	25.3	3.8	0.0	24.0	0.0	0.0
cy.20	46.1	24.2	7.7	34.9	16.3	7.0
cy.21	38.7	2.7	0.0	37.0	0.0	0.0
cy.22	31.3	3.0	0.0	21.0	0.0	0.0
cy.23	31.4	5.9	0.0	16.7	0.0	0.0
cy.24	35.4	8.3	0.0	⋯	⋯	⋯

Note. The values show the percentages of the observed sunspot number, which are out of the 1σ, 2σ, and 3σ of the predicted variations.

Download table as:  ASCIITypeset image

Figure 6 shows the comparisons between the predicted and the observed sunspot number and the latitudinal location of sunspot groups (discussed in the next section). We take the following cases as examples: (1) cycle 20 (first row) and cycle 23 (second row) predicted at 4 years into cycles and showing smaller goodness-of-fits than HWR94, and (2) cycle 12 (third row) and cycle 14 (fourth row) predicted at 8 years into the cycles and showing larger goodness-of-fits than HWR94. The first column is the smoothed monthly total sunspot number. Although the mean values of our predictions still always show the deviations from the observations, the observations are usually within our predicted 2σ uncertainty ranges, denoted by the shaded regions. This also demonstrates the importance of giving the uncertainty range. The second column is the time evolution of the monthly sunspot number. The random realizations look very similar to the observed ones.

Zoom In Zoom Out Reset image size

3.2. Prediction of the BMR Emergence

Once we have the predicted sunspot number, the next step is to predict the emergence properties of the sunspot groups. Jiang et al. (2011) give the dependence of the statistical properties of sunspot group emergence, including the latitude, longitude, area, and tilt angle of sunspot groups in the form of the BMRs on the cycle phase and strength. Those properties are used here. The following list summarizes the statistical properties aiming to convert the sunspot number time series into the time series of sunspot group emergence.

• 1.  
  The number of BMRs emerging per month was taken to be equal to ${F}_{\mathrm{BMR}}(t)=0.24{F}_{\mathrm{mn}}(t)$. The number of the daily emergence of the BMR is randomly realized, satisfying the total number F_BMR(t).
• 2.  
  The mean latitudinal distribution of BMR emergence ${\lambda }_{n}$ depends on the cycle strength S_n. It obeys the form ${\lambda }_{n}(x)=(26.4-34.2x+16.1{x}^{2})({\overline{\lambda }}_{n}/\langle {\lambda }_{n}{\rangle }_{12-20})$, where ${\overline{\lambda }}_{n}=(12.2+0.015{S}_{n})/14.6$ and x is the fraction of a solar cycle period.
• 3.  
  The width of the latitude distribution obeys a Gaussian profile with a half width of σ, which is equal to $(0.14+1.05x-0.78{x}^{2}){\lambda }_{n}(x)$. We exclude points deviating from the mean by more than 2.2σ for the equatorward side.
• 4.  
  The BMR emergence has the symmetric distribution in north and south hemispheres.
• 5.  
  The BMR emergence has the random longitudinal distribution.
• 6.  
  The number density function of sunspot group areas obey the following distributions:

  $$\begin{eqnarray}&&\begin{array}{l}{n}_{\mathrm{BMR}}({A}_{S})\\ =\left\{\begin{array}{ll}0.3{A}_{S}^{-1.1} & \mathrm{for}\,{A}_{S}\leqslant 60\mu \mathrm{Hem}\\ 0.003\exp [-\displaystyle \frac{1}{2\mathrm{ln}3}{(\mathrm{ln}{A}_{S}-\mathrm{ln}45)}^{2}] & \mathrm{for}\,{A}_{S}\gt 60\mu \mathrm{Hem}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 19 }$$
• 7.  
  The average areas of sunspot groups are cycle phase (x) dependent, which obeys $\overline{{A}_{S}}(x)=115\,+\,396x\,-426{x}^{2}$.
• 8.  
  The umbra area of each BMR is ${A}_{U}={A}_{S}/5$. The total area of BMR $A={A}_{S}+{A}_{P}$, where A_P is the plage area and is calculated by ${A}_{P}=414+21{A}_{S}-0.0036{A}_{S}$ (Chapman et al. 1997).
• 9.  
  The mean tilt angle, $\overline{{\alpha }_{n}}(\lambda )$, of the emerging BMRs for the cycle n is assumed to follow Joy's law in the form $\overline{{\alpha }_{n}}={T}_{n}\sqrt{| \lambda | }$, where ${T}_{n}=1.72-0.0022{S}_{n}$.
• 10.  
  The averaged value of the scatter of the tilt angle, ${\rm{\Delta }}\alpha$, is zero and its standard deviation is ${\sigma }_{\alpha }$. We use the empirical relation between ${\sigma }_{\alpha }$ and A_U, which is that ${\sigma }_{\alpha }=-11\mathrm{log}{A}_{U}+35$, to derive ${\sigma }_{\alpha }$ (Jiang et al. 2014a).
• 11.  
  A factor of 0.7 is multiplied to the resulting tilt angle ${\alpha }_{n}=\overline{{\alpha }_{n}}+{\sigma }_{\alpha }$ to account for the effect of latitudinal inflow toward BMRs (Jiang et al. 2010; Martin-Belda & Cameron 2017).

The third column of Figure 6 shows examples of the time evolution of the sunspot latitudinal distribution, i.e., butterfly diagrams. The predicted diagram is consistent with the observed one.

Zoom In Zoom Out Reset image size

The BMR emergence provides the source of the large-scale magnetic field over the solar surface. The evolution of the large-scale magnetic field can be described by the SFT model, which will be presented in the following section.

4.1. SFT Model

The relevant equation to describe the evolution of the large-scale magnetic flux distribution at the solar surface $B(\lambda ,\phi ,t)$ as a combined result of the emergence of BMRs, a random walk due to supergranular flows, and the transport by large-scale surface flows is

$$\begin{eqnarray}&&\displaystyle \frac{\partial B}{\partial t}=-{\rm{\Omega }}(\lambda )\displaystyle \frac{\partial B}{\partial \phi }-\displaystyle \frac{1}{{R}_{\odot }\cos \lambda }\displaystyle \frac{\partial }{\partial \lambda }[\upsilon (\lambda )B\cos \lambda ]\\ &&\quad +\,\eta \left[\displaystyle \frac{1}{{R}_{\odot }^{2}\cos \lambda }\displaystyle \frac{\partial }{\partial \lambda }\left(\cos \lambda \displaystyle \frac{\partial B}{\partial \lambda }\right)+\displaystyle \frac{1}{{R}_{\odot }^{2}{\cos }^{2}\lambda }\displaystyle \frac{{\partial }^{2}B}{\partial {\phi }^{2}}\right]\\ &&\quad +\,S(\lambda ,\phi ,t),\end{eqnarray} \tag{ 20 }$$

where λ and ϕ are the heliographic latitude and longitude, respectively. The latitudinal differential rotation of Sun, ${\rm{\Omega }}(\lambda )$, is taken from Snodgrass (1983) and the poleward meridional flow, $\upsilon (\lambda )$, is given by van Ballegooijen et al. (1998) and is consistent with the measurement by Hathaway & Rightmire (2011). The magnetic diffusivity, η, is 250 km² s⁻¹, which is within the range from the observational studies summarized by Schrijver et al. (1996). All of the three transport parameters are time independent. $S(\lambda ,\phi ,t)$ is the time-dependent flux source term, which will be obtained based on the predicted sunspot emergence using the empirical method in Section 3. The initial magnetic field configuration due to the BMR emergence based on its area, location, and tilt is given by Baumann et al. (2004). The corresponding magnetic flux is determined by a single parameter, Bmax (=592 G), which is calibrated by the total unsigned surface flux obtained from SOHO/MDI polar-field-corrected synoptic maps (Sun et al. 2011) after rebinning the spatial resolution of the simulation (1° in both latitude and longitude).

For the initial condition of each SFT simulation, we use the SOHO/MDI polar-field-corrected synoptic maps (available from 1996 June until 2010 May) and the radial synoptic maps with 3600 points in Carrington longitude and 1440 points equally spaced in sine latitude from HMI (2010 May until the present). The abnormal points that correspond to unobserved values or the absolute values over 50 G above $\pm 60^\circ$ latitudes are filled by the averaged value of neighboring normal 3 points. Both the MDI and the HMI data are smoothed with the width of 7. Both data were reduced to a resolution of 1° in latitude and longitude by the IDL CONGRID function and then were converted to the equal latitudes. The rightmost panel of Figure 8 shows that the cross-calibration of HMI and MDI based on the axial dipole moment during the overlap time period (Carrington Rotations, CRs2097-2104) is 1.3. This is consistent with Liu et al. (2012) who compared line-of-sight magnetograms taken by MDI and HMI. We multiply by a factor 1.3 when the HMI synoptic map is used as the initial condition. The low-resolution data of synoptic maps were used by Cameron & Schüssler (2016) and Jiang & Cao (2018).

Zoom In Zoom Out Reset image size

We used the code originally developed by Baumann et al. (2004) to do the numerical calculations. The magnetic field is expressed in terms of spherical harmonics up to l = 63. A fourth-order Runge–Kutta method is used for time stepping with an interval of one day.

4.2. Prediction of the Large-scale Magnetic Field Evolution over Surface

With 50 sets of random realizations of the magnetic flux emergence and the synoptic magnetograms as the initial condition, we may derive the possible large-scale magnetic field evolution over the surface using the SFT model. Since here our interest mainly concentrates on long-term predictions, we mainly show the results about the time evolution of the axial dipole moment D(t), which is defined as

$$\begin{eqnarray}&&D(t)=\displaystyle \frac{3}{2}{\int }_{0}^{180}\langle B\rangle (\theta ,t)\cos \,\theta \sin \,\theta d\theta .\end{eqnarray} \tag{ 21 }$$

The uncertainties of the results originate from two ingredients. One is the uncertainty due to the scatter in the properties of the BMRs' emergence, which are assumed to be reasonably represented in Sections 2 and 3. The other is the uncertainty due to measurement error in the magnetograms used as initial conditions. The imperfect measurements mainly result from the instrumental characteristics, e.g., spatial resolution, scattered light, and filter characteristics, and from the inherent complexity of the solar magnetic field, e.g., the saturated factors for different spectral lines (Wang & Sheeley 1995). The properties of the magnetic field at lower latitudes have larger effects on the results (Jiang et al. 2014a). Since some imperfections of the measurements are inevitable, they are potential uncertainties to the predictive skill of the model. Here we only consider the averaged radial surface field over the surface, which deviates from zero due to measurement error. We use the same method described in Cameron et al. (2016) to estimate the error. The total error due to the magnetogram measurement and the random flux source emergence is determined by adding them quadratically.

We take the timing of 4 and 8 years into cycles 23 and 24 as examples. Figure 7 shows the predicted time evolutions of the axial dipole moments for cycle 23 and cycle 24 at 4 and 8 years into the cycles. Synoptic magnetograms CR1961 and CR2013 from MDI and CR2130 and CR2184 from HMI are used as the initial magnetic field, respectively. Solid green lines show the average of 50 SFT simulations with random sources. Dark and light red shading indicate the total σ and 2σ uncertainties, which correspond to the probabilities of 68% and 95.4% for the observed axial dipole moment to be within the error ranges, respectively. The dashed green lines give the 2σ range for the intrinsic solar contribution (source scatter). The errors from MDI measurements are significant and are larger than that from HMI measurements. This also demonstrates the big effects of observed magnetograms on the prediction.

Figure 7 shows the following results. First, the observed values are always within the 2σ uncertainty, although statistically it is possible that the observation is out of the 2σ range with the probability of 4.6%. This demonstrates the effectiveness of prediction method. Second, the error ranges increase with time. The error ranges of the predicted axial dipole moment at the end of the cycles are smaller when the prediction is made later. This is due to the accumulative effect of the random source emergence. Third, the absolute values of the predicted mean axial dipole moments decrease when the prediction was made at a later time. Cycle 23 had a deep minimum with very low axial dipole moment. Jiang et al. (2015) found that the weaknesses of the cycle 23 minimum are mainly caused by a number of bigger bipolar regions emerging at low latitudes with a wrong (i.e., opposite to the majority for this cycle) orientation of their magnetic polarities in the north–south direction, which impaired the growth of the polar field. The abnormal emergences are supposed to be emerged during cycle 24 as well. Hence the mean value at cycle minimum decreases with time in both cycle 23 and cycle 24. The corresponding predictions of the following cycle strength decrease with time as well. But it is possible that the predicted mean values of the axial dipole moment increase with time if the peculiar BMR emergences (large areas and tilts at low latitudes) with positive contributions to the axial dipole moment surpass the peculiar BMR emergences with negative contributions to the axial dipole moment. In such cases, a stronger subsequent cycle will be proceed.

The prediction of the subsequent cycle depends on the correlation between the dipole moment at cycle minimum and the following cycle strength. The physical base has been clarified in Section 1. This correlation is a key ingredient to affect the prediction accuracy of the subsequent cycle. A long-term homogeneous dipole moment data set, which is directly calculated based on the observed synoptic magnetograms, is expected. The earliest available synoptic magnetograms are from Mount Wilson Observatory (MWO) from 1974 July until 2012 December.⁷ Another long-term data set that is widely used in solar cycle prediction is the magnetic measurement from Wilcox Solar Observatory (WSO, continuously available from 1976 April until present⁸ ). We use the axial dipole moment calculated based on MDI synoptic maps to cross-calibrate the other three axial dipole moments calculated based on MWO, WSO, and HMI synoptic maps.

Figure 8 shows the cross-calibrations of the axial dipole moment from MDI (${D}_{\mathrm{MDI}}$) with that from MWO (${D}_{\mathrm{MWO}}$, left panel), from WSO (${D}_{\mathrm{WSO}}$, middle panel) and from HMI (${D}_{\mathrm{HMI}}$, right panel). The linear fits give their relations ${D}_{\mathrm{MDI}}=0.17+3.98{D}_{\mathrm{MWO}}$, ${D}_{\mathrm{MDI}}=0.14+5.42{D}_{\mathrm{WSO}}$, and ${D}_{\mathrm{MDI}}=1.3{D}_{\mathrm{HMI}}$. The correlation coefficients are 0.94, 0.96, and 0.89, respectively.

Figure 9 is the time evolution of the calibrated axial dipole moment after 13 months of running averages in solid curves overplotted with sunspot number shown by the dashed curve from 1974 onwards. We derive the time series of the averaged axial dipole moment $\overline{D}(t)$ in the following way. If there is more than one observation, we determine the average over the different values. During the early time, we just take the values from MWO observations since only the MWO observation is available. The 13-month running average of $\overline{D}(t)$ is regarded as the homogeneous axial dipole moment data set to derive the correlation between the dipole moment at cycle minimum ${D}_{\min }^{n}$ and the following cycle strength ${S}_{n+1}$. To further reduce the random noise, we use the averaged value of 7 CRs around each minimum to get ${D}_{\min }^{n}$. The cycle strength ${S}_{n+1}$ is the average of 7 months around each maximum of the smoothed sunspot data. The correlation between the homogeneous dipole moment at cycle minima ${D}_{\min }^{n}$ and the following cycle strength ${S}_{n+1}$ is shown in Figure 10, which indicates

$$\begin{eqnarray}&&{S}_{n+1}=58.7\,\ast \,{D}_{\min }^{n}.\end{eqnarray} \tag{ 22 }$$

Although there are only four points, the correlation coefficient is r = 0.99 with the corresponding confidence level p = 0.045. The sudden jump of the WSO observations during 2016–2017 in Figure 9 is due to the reduced WSO polarization sensitivity.⁹ This time period has no effect on the correlation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Prediction of the Subsequent Cycle at Different Phases of a Cycle

We can predict the possible sunspot emergence of an ongoing cycle by random realizations of the features of sunspot emergence given in Section 2 a few years after the start of a cycle. With the SFT simulations and the observed synoptic magnetograms as the initial state of the magnetic field, the possible large-scale field evolution over surface can be obtained. Based on the correlation between the dipole moment at cycle minimum and the following cycle strength constrained in Section 5 combined with the features of sunspot emergence given in Section 2, the predictability of the shape of the subsequent cycle can be given.

Solar cycles also vary in length. Strong cycles tend to be shorter and vice versa, but with a large scatter (Solanki et al. 2002; Du 2006; Vaquero & Trigo 2008). Some of this is due to the overlapping of sunspots from adjacent cycles. In this paper, we ignore the cycle length variation and assume that all the cycles have an 11-year cycle period. If we know the maximum of a cycle S_n, we may derive the profile of the cycle based on the simplified formula of HWR94, which is just a function of a relevant to cycle amplitude S_n. The parameter t₀ is not required. The correction function $\overline{\bigtriangleup r(t)}$ is still required for the improvement of the later phase. We only predict the systematic part of the smoothed sunspot number ${\overline{{F}_{\mathrm{sm}}}}^{S}(t)$ for the prediction of the subsequent cycle, since the standard deviation of the prediction is usually larger than the intrinsic short-term variability. The uncertainty results from the uncertainty of the axial dipole moment at the end of the cycle.

The formula to describe the subsequent cycle is

$$\begin{eqnarray}&&{\overline{{F}_{\mathrm{sm}}}}^{S}(t)={f}_{2}(t)+\overline{\bigtriangleup r(t)}{f}_{2}(t),\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&{f}_{2}(t)=\displaystyle \frac{{{at}}^{3}}{\exp [{t}^{2}/{b}^{2}]-c},\end{eqnarray} \tag{ 24 }$$

and $\overline{\bigtriangleup r(t)}$ corresponds to Equation (11). The fits of cycles 12–24 using Equation (9) indicate that the relationship between S_n and a is

$$\begin{eqnarray}&&{S}_{n}=9072.8{a}^{0.706}.\end{eqnarray} \tag{ 25 }$$

Parameters b and c are the same as those in Equation (1). The deviation of the maximum sunspot number given by Equation (24) from S_n is within 5.

Figure 11 shows some examples illustrating the predictability of the subsequent cycle, starting from different phases of a cycle. The timings for the predictions are the same as those in Figure 7. From Figure 7, we can derive ${\overline{D}}_{n}+2{\sigma }_{n}^{D}$, ${\overline{D}}_{n}+{\sigma }_{n}^{D}$, ${\overline{D}}_{n}$, ${\overline{D}}_{n}-{\sigma }_{n}^{D}$, and ${\overline{D}}_{n}-2{\sigma }_{n}^{D}$ at the cycle n minimum. We derive the corresponding amplitude of the subsequent cycle strength, ${\overline{S}}_{n+1}+2{\sigma }_{n+1}^{S}$, ${\overline{S}}_{n+1}+{\sigma }_{n+1}^{S}$, ${\overline{S}}_{n+1}$, ${\overline{S}}_{n+1}-{\sigma }_{n+1}^{S}$, and ${\overline{S}}_{n+1}-2{\sigma }_{n+1}^{S}$ using Equation (22). Still, the ${\sigma }_{n+1}^{S}$ and $2{\sigma }_{n+1}^{S}$ ranges correspond to the probabilities of 68% and 95.4% for the actual sunspot number to be within the error ranges. Based on Equations (23)–(25), profiles of the smoothed sunspot number in the subsequent cycle $n+1$ can be obtained. We see that the observed sunspot number in cycle 24 is within the σ to 2σ range of the predicted result, whatever the prediction is at the early time, i.e., the fourth year or the later time, i.e., the eighth year. The error range decreases with time. The deviation between the observation and prediction during the decline phase is expected since the prediction given here is only the systematic part. As shown in Table 1, the systematic part sometimes shows a deviation from the observed sunspot number evolution (but the observations still are within the given error range). The updated prediction of the cycle based on a later magnetogram and the random flux source is presented in the following subsection.

Zoom In Zoom Out Reset image size

6.2. Updated Predictability of Cycle 25

Here we take the synoptic magnetogram CR2198 (2017 December 2nd–31st) as the initial condition to investigate the large-scale field evolution over the solar surface during the rest of cycle 24 and the possible profiles of cycle 25. A possible (one from our Monte-Carlo ensemble) monthly sunspot emergence pattern during the rest of the ongoing cycle 24 is shown in Figure 12. The systematic part of the smoothed sunspot number is shown in the solid green curve. It is higher than the current observation. But the deviations are within the 2σ range of the fluctuation. The thin green curve is one realization of the smoothed sunspot number. The corresponding time evolution of the latitudinal location of the sunspot groups (observed in black and predicted in red) is presented in the upper left panel of Figure 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The upper right panel of Figure 13 shows the time evolution of the axial dipole moment from MDI and HMI observations in blue and in black before the end of 2017 and the predicted one after 2017. The mean value will increase from 1.66 G at 2018.0 to 2.12 G with a 2σ range of 0.55 G at 2020. We also calculated the averaged polar field over $\pm 60^\circ$ to $\pm 75^\circ$ latitudes to understand the different evolutions of the two hemispheres. The northern polar field is the dashed curve and the southern polar field is the solid curve. The green curves show the expected values and the dark and light red shades show the 1σ and 2σ range of the predicted polar field, respectively. During the first year, the error range is very small. The polar field is determined by the initial conditions. The northern and the southern polar fields are almost in balance with values of 3.17 and −2.86 G. The southern polar field will keep almost flat with values varying from −2.94 G, and the northern polar field will have a large increase to 4.56 G at the end of 2018. The averaged northern polar fields keep increasing and the averaged southern polar fields keep stable until the end of the cycle. The mean values by then are −2.89 G and 5.40 G with a 2σ range of 1.23 G, respectively. The reason for the increasing northern polar field and the stable southern polar field can be explained by the lower right panel of Figure 13, which is the longitudinal averaged surface field evolution combined with the HMI observations (before the vertical line) and the simulation (after the vertical line) corresponding to the random realization in upper left panel. There are strong poleward positive plumes starting from second half of 2017. The positive plume in the northern hemisphere further increased the positive polar field. The positive plume in the southern hemisphere prevented the increase of southern negative polar field. Hence it keeps almost stable. The two notable plumes mainly result from two big ARs, i.e., AR12674 and AR12673 (Yang et al. 2017; Yan et al. 2018). They occurred on the solar disc with large tilts around 2017 September 5th located at northern and southern hemispheres, respectively. According to Jiang et al. (2014a, 2015), the bigger bipolar regions emerging at low latitudes have significant effects on large-scale field evolution. Their effects on the solar cycle will be studied in detail in a separated study.

The possible evolution of the smoothed sunspot number based on the axial dipole moment at the end of cycle 24 is given in Figure 12. The expected maximum amplitude of cycle 25 is 125, which is about 10% higher than current cycle 24. The 2σ range is 32, which means the possibility that the amplitude of cycle 25 above 93 is 95.4%. Cycle 25 most probably is a normal cycle, rather than the Maunder minimum period. The possible profiles of cycle 25 denoted by the shaded region obey the empirical relations given in Section 2.

In order to avoid effects due to the unexpected problems from synoptic magnetograms, we also repeated the predictions based on synoptic magnetograms CR2196 and CR2197. The results are similar.

In this paper, we have developed a scheme to investigate the predictability of the solar cycle over one cycle. The scheme includes three steps. First, empirical properties of the solar cycle are used to predict the possible sunspot emergence for an ongoing cycle. Then the SFT model is adopted to predict the possible large-scale field evolution over the surface, including the polar field at the end of the cycle. Finally, the correlation between the polar field and the subsequent cycle strength and empirical properties of the sunspot emergence are applied to get the possible profiles of the subsequent cycle. The scheme is verified by past cycles and is applied to predict the possible profiles of cycle 25. The results show the cycle 25 strength of 125 ± 32 (2σ uncertainty range), which is about 10% stronger than cycle 24 based on the mean value.

Comparing with the existing methods of the solar cycle prediction, the main progress of the current scheme is as follows. First, the prediction scope is extended over one cycle. Second, the profiles of solar cycles during the whole prediction time period, rather than only the cycle amplitude, can be given. Third, the uncertainty due to the randomness of the sunspot emergence and the data assimilated to the model during the whole prediction time period is given. Fourth, not only the monthly and smoothed sunspot number but also the emergence properties, including the size, location, and tilt of sunspot groups, can be given. Fifth, although the BL dynamo model is not directly used in the model, the empirical properties we used have the solid dynamo origin. These five properties distinguish the scheme from all of the current existing attempts of the solar cycle predictions.

For one source of prediction uncertainties, the imperfect measurement of the initial magnetogram, we only estimate the error from the net flux density. The error due to other possible problems, like the center-to-limb correction of the magnetic saturation, is not dealt with in this paper. Potential problems of the initial magnetograms can degrade the predictive ability. We have explicitly assumed that the scatter in the tilt angle is random in nature, consistent with the idea that it is due to turbulent buffeting by the convective motions. If there is a significant deterministic chaotic component in the dynamo process, which we currently do not see evidence for, the model can be improved. Conversely if there is some chaotic or random process operating on long timescales and which exceeds the randomness introduced by the scatter in the tilt angles, then our prediction can fail. Furthermore, we assumed that all cycles have an 11-year cycle period and ignored the cycle overlap. So there remains room for improving the scheme, although these improvements are likely to be small compared to the uncertainties due to the tilt-angle scatter that we explicitly dealt with in this paper.

We are grateful to Robert Cameron, Sami Solanki, and Manfred Schüssler for stimulating discussions. Robert Cameron also helped to improve the language. SOHO is a project of international cooperation between ESA and NASA. The SDO/HMI data are courtesy of NASA and the SDO/HMI team. The sunspot records are courtesy of WDC-SILSO, Royal Observatory of Belgium, Brussels. The National Solar Observatory (NSO)/Kitt Peak data used here are obtained cooperatively by NSF-NOAO, NASA/Goddard Space Flight Center, and the NOAA Space Environment Laboratory. NSO/SOLIS data were courtesy of NISP/NSO/AURA/NSF. We acknowledge the support by the National Science Foundation of China (grant Nos. 11522325, 11573038, and 41431071) and by the Fundamental Research Funds for the Central Universities of China.