Modelling the evolution of the Sun's open and total magnetic flux

Solar activity in all its varied manifestations is driven by the magnetic field. Particularly important for many purposes are two global quantities, the Sun's total and open magnetic flux, which can be computed from sunspot number records using models. Such sunspot-driven models, however, do not take into account the presence of magnetic flux during grand minima, such as the Maunder minimum. Here we present a major update of a widely used simple model, which now takes into account the observation that the distribution of all magnetic features on the Sun follows a single power law. The exponent of the power law changes over the solar cycle. This allows for the emergence of small-scale magnetic flux even when no sunspots are present for multiple decades and leads to non-zero total and open magnetic flux also in the deepest grand minima, such as the Maunder minimum, thus overcoming a major shortcoming of the earlier models. The results of the updated model compare well with the available observations and reconstructions of the solar total and open magnetic flux. This opens up the possibility of improved reconstructions of sunspot number from time series of cosmogenic isotope production rate.


Introduction
The precise history of solar activity and its underlying magnetic field is of interest for a number of reasons. Firstly, records of solar activity and of the magnetic field pose an important constraint on models for the enhancement and the evolution of magnetic flux (mainly dynamo models). Secondly, such records are important for understanding the history of the Sun's influence on the Earth (either through changes in its irradiance, or through space weather effects). Thirdly, a long record of solar activity is needed to understand how the Sun compares with other sunlike stars in its level of activity and variability (e.g., Radick et al. 2018;Reinhold et al. 2020).
Solar activity, in all its diverse manifestations, is driven by its magnetic field, so that knowledge of the history of solar activity implies knowledge of its magnetic field. Two widely used quantities describing the global magnetic field of the Sun are the global open and unsigned total magnetic flux. They are, for example, used in heliospheric physics, for the reconstruction of solar irradiance, or as measures of solar activity when comparing with other stars. These quantities, being global in nature, can be reconstructed from more indirect proxies of solar activity and magnetism, such as the sunspot number and concentrations of the cosmogenic isotopes 14 C or 10 Be in terrestrial archives.
A first model to compute the solar open flux from the sunspot number was developed by Solanki et al. (2000), based on a simple differential equation describing the evolution of the open flux, φ open . In spite of its simplicity, it successfully reproduced the empirically-reconstructed evolution of the open flux by Lockwood et al. (1999) and the 10 Be concentration in ice cores (Beer et al. 1990). This simple model was extended by Solanki et al. (2002) to cover also the total unsigned magnetic flux, φ total , but now requiring the solution of a set of coupled differential equations to describe the evolution of ephemeral re-gions (ERs) besides that of active regions (ARs) and of the open flux. All three components of magnetic flux contribute to the evolution of the total magnetic flux. ARs are the large bipolar structures that harbour sunspots at least part of the time, whereas ERs are smaller bipolar regions without sunspots.
The ability of this model to reproduce concentrations of cosmogenic isotopes turned out to be particularly useful (e.g., Usoskin et al. 2002). Although far more sophisticated models are in the meantime available to compute not just global magnetic quantities, but also the underlying spatial distribution of the magnetic flux and the detailed input from individual emerging ARs etc., the very simplicity of this set of models allowed them to be inverted (e.g., Lockwood 2003), so that, e.g., sunspot number could be reconstructed from measured concentrations of cosmogenic isotopes (e.g., Usoskin et al. 2003Usoskin et al. , 2004Solanki et al. 2004;Usoskin et al. 2016;Wu et al. 2018b). The model of Solanki et al. (2002) was further extended and combined with the successful SATIRE model (Spectral And Total Irradiance REconstruction; Fligge et al. 2000;Krivova et al. 2003Krivova et al. , 2011 to compute total solar irradiance over the last 400 years (Krivova et al. 2007(Krivova et al. , 2010.  have further refined the model by distinguishing between the short-lived and long-lived components of the open flux (Ikhsanov & Ivanov 1999;Cranmer 2002;Crooker et al. 2002;see Vieira & Solanki 2010 for details), which led to an improved reconstruction of the open flux that displayed a better agreement with observations. This model, with some tuning, has been the basis for further reconstructions of solar spectral irradiance over the telescope era (Krivova et al. 2010), as well as sunspot number and TSI over the Holocene (Vieira et al. 2011;Wu et al. 2018b,a).
One shortcoming of earlier versions of the model discussed above is that the open flux during a grand minimum, such as the A&A proofs: manuscript no. AA_2021_40504 Maunder minimum, i.e. during a long period essentially without sunspots, invariably drops to zero. This is because in this model, the emergence rate of the magnetic field on the solar surface is linearly linked to the sunspot number, so that by design during a grand minimum no magnetic flux is allowed to emerge. This leads to a zero open and total flux during the grand minima. It has been shown, however, that signals of solar activity and variability were also present during the Maunder minimum (Beer et al. 1998;Fligge et al. 1999;Usoskin et al. 2001;Miyahara et al. 2004;Riley et al. 2015). This was also confirmed by modelling (Owens et al. 2012) and points to a need for an improvement of the global total and open magnetic flux model. Furthermore, more recent solar observations provided new insights into the sources, emergence and evolution of the solar magnetic flux. Thus, Thornton & Parnell (2011) have combined observations from various sources and found that the emergence rate of bipolar magnetic regions with fluxes between 10 16 Mx and 10 23 Mx follows a single power law with a slope of −2.69.
Here we present a new, strongly revised version of the VS2010 model that builds on these recent solar observational results, replacing the direct proportionality of ERs and ARs by a more up-to-date approach. It does keep the original differential equations, however, so that it is not a completely independent model. As a natural outcome of the model, ERs keep emerging even during a grand minimum when there are no sunspots for multiple decades. This means that neither the open nor the total magnetic flux drop to zero at any time.
The paper is structured as follows. The data used to constrain and test the new model are briefly introduced in Sect. 2. We describe our model and highlight the changes relative to the older version of the model in Sect. 3. The results of the model are presented in Sect. 4, while we summarise and discuss our findings in Sect. 5, where we also provide an outlook on future applications of the new model.

Data
The model to be detailed in Sect. 3 starts from a sunspot number time series and computes the total and open magnetic fluxes of the Sun therefrom. We therefore require a sunspot series as input to the model. To constrain the free parameters of the model, we compare its output to observations and independent data-based reconstructions of the total and open magnetic fluxes. Finally, to test the output of the model, we consider further independent time series of the reconstructed open magnetic flux.
As input to the model we use the following sunspot number data sets: (1) the international sunspot number v2.0, referred to hereafter as ISN2.0 (Clette & Lefèvre 2016), and (2) the group sunspot number, or GSN in short (Hoyt & Schatten 1998). The ISN2.0 data set was extended back to 1643 by adding the sunspot data during the Maunder minimum by Vaquero et al. (2015) scaled up by the factor 1.67 to match the ISN2.0 definition. Also in the GSN record, the values before 1710 were replaced by the data from Vaquero et al. (2015), without any scaling.
To constrain the output of the model we make use of observations of the total magnetic flux (see Arge et al. 2002;Wang et al. 2005;Wenzler et al. 2006) derived from synoptic charts produced by the three solar observatories with the longest running regular magnetographic measurements: Wilcox Solar Observatory (WSO), Mount Wilson solar Observatory (MWO), and National Solar Observatory at Kitt Peak (NSO/KP). These data sets have already been used to constrain earlier versions of this model and we have used the same versions as employed by Krivova et al. (2007Krivova et al. ( , 2010 and Wu et al. (2018a). To constrain the free parameters of the model, we consider the average over at least two (MWO andWSO during 2002-2009) or all three  records for each Carrington rotation.
Furthermore, to better constrain the free parameters of the model, we also use the empirical reconstruction of the open magnetic flux from the geomagnetic aa-index covering the period from 1845 to 2010 (Lockwood et al. 2014).
Finally, the quality of the computed open magnetic flux is tested by comparing it with two other independent data sets (without changing the free parameters of the model): (1) a compilation of spacecraft-based in-situ measurements by Owens et al. (2017) since 1998, and (2) a reconstruction by Wu et al. (2018b) from decadal INTCAL13 14 C data covering the Holocene prior to 1900 (Reimer et al. 2013).

Magnetic flux emergence and evolution
Following the approach by Solanki et al. (2000Solanki et al. ( , 2002 and , we describe the evolution of the solar total and open magnetic flux by a system of ordinary differential equations. However, instead of distinguishing between active regions (AR) and ephemeral regions (ER) as the sources of fresh magnetic flux at the solar surface, we distinguish between ARs and what we call Small-Scale-Emergences (SSEs), i.e. all emergences with fluxes smaller than those in active regions. These therefore combine the flux emerging in the form of ephemeral regions and smaller magnetic bipoles all the way down to internetwork fields. The equations describing the evolution of the different (globally averaged) components of the magnetic flux are: Here In the original model, the emergence rate of ARs at a given time t, ε AR (t), was linked linearly to the sunspot number, SN, at that time: where ε max,21 AR and SN max,21 are the three-month averaged emergence rate and SN value observed during the maximum of cycle 21 (taken from Schrijver & Harvey 1994), respectively. Because at the time that the model was originally developed, large-scale studies of magnetic flux emergence and evolution could not resolve internetwork fields, the earlier model was restricted to ERs as the only magnetic bipoles smaller than ARs. The emergence rate of the ERs, ε ER,n , of the cycle n was not well known and was assumed to be a sinusoidal function, g n , which was stretched such that the length of the ER cycle was longer than the respective sunspot cycle. The amplitude of the ER cycle was simply taken to be proportional to the maximum value of the emergence rate of the ARs of the cycle, so that the emergence rate of ERs over a cycle had the form: where X is an amplitude factor (a free parameter of the model; same for all cycles). Importantly, Eq. 8 implies a linear relationship between the emergence rate ε ER,n and the sunspot number. This, along with the linear relationship between ε AR,n and sunspot number (Eq. 7) results in an absence of flux emergence during extended periods of spotless days.  (9) following Thornton & Parnell 2011). φ 0 , φ ER , φ AR and φ limit represent the limit below which the local dynamo flux dominates, the minimum ephemeral region flux, the minimum active region flux, and the upper limit of the active region flux, respectively. The horizontal arrows mark the flux ranges corresponding to the internetwork (IN), ephemeral regions (ERs), active regions (ARs), and Small-Scale Emergences (SSEs). The SSE range includes IN fields and ERs. The slope m of the distribution was derived by Thornton & Parnell (2011) by fitting various observations at different activity levels. Slopes m 1 and m 2 represent the corresponding distributions at maximum and minimum of solar activity levels for cycle 21, respectively (see main text for details).
To overcome this shortcoming, we incorporate more recent solar observations by Thornton & Parnell (2011) into the model. Using high-resolution Hinode, Solar Optical Telescope/Narrowband Filter Imager (SOT/NFI) observations and combining them with earlier published data, they found that the emergence rate of the magnetic flux on the Sun follows a power-law distribution: where φ 0 = 10 16 Mx (the smallest flux per feature that they include in their histograms), n 0 = 3.14 × 10 −14 cm −2 day −1 and m = −2.69 (see the illustration in Fig. 1). We note that Thornton & Parnell (2011) have summarised the results from multiple studies with a wide range of solar activity levels and observing conditions. In earlier studies, Harvey & Zwaan (1993) and Harvey (1993) found that the emergence rate of ARs varied significantly more between solar activity minimum and maximum than that of ERs. Whereas roughly 8.3 times more ARs emerged during the maximum of cycle 21 than during the minimum (the factor generally grows with the size of the regions, but was on average about 8.3 for all the ARs they studied), this ratio was roughly two for ERs. The number of the smallest magnetic features, forming the internetwork magnetic fields and having fluxes of 10 16 -10 17 Mx, appears to be nearly invariable over an activity cycle (Buehler et al. 2013;Lites et al. 2014). These features differ from the larger ones in that they are mainly brought about by a small-scale turbulent dynamo (Vögler & Schüssler 2007;Rempel 2014) that produces the same amount of magnetic flux nearly independently of largescale activity.
To satisfy these observational constraints, on the one hand, we keep the number of the smallest magnetic features considered here (with a flux per feature of 10 16 Mx) fixed at all times. On the other hand, we allow the exponent m to vary (cf. Parnell et al. 2009;Schrijver & Harvey 1989) around the empirical value m = −2.69 found by Thornton & Parnell (2011) The slopes m 1 and m 2 describe the distributions of emergence rates during periods when the observed sunspot numbers are SN 1 and SN 2 (with SN 1 > SN 2 ). In our model, the slope m follows the SN, m(SN), according to the non-linear relationship: where α is a free parameter, fixed by comparing the output of the model to observations and independent reconstructions (see Sect. 3.2). Now, the emergence rate of magnetic flux in ARs and SSEs at any given time can be calculated as: and Here φ AR = 4 × 10 20 Mx denotes the magnetic flux of the smallest bipolar regions hosting sunspots, i.e. the smallest active regions (e.g., Zwaan 1978;Schrijver & Zwaan 2000; van Driel-Gesztelyi & Green 2015) and φ limit is the flux of the largest considered ARs. Since such regions are extremely rare, the exact value of φ limit is not important. Following Parnell et al. (2009) and Thornton & Parnell (2011), we take φ limit = 10 23 Mx, which is somewhat larger than the maximum flux (3 × 10 22 Mx) for ARs listed by Schrijver & Zwaan (2000) and van Driel-Gesztelyi & Green (2015). Tests have shown that also the exact value of φ AR adopted here plays only a minor role for the end result in the sense that although the free parameters To estimate how the slope m of the distribution given by Eq. (9) changes with activity (Eq. (10)) we rely on the observations by Harvey & Zwaan (1993) and Harvey (1993) for cycle 21. They found that the number of emerging ARs in cycle 21 varied between the activity maximum and minimum by a factor of roughly 8.3. The monthly-smoothed sunspot numbers corresponding to these periods (1979 -1982 for the maximum, as well as 1975 -1976 and 1985 -1986 for the preceding and following minima) are then SN 1 = 217 and SN 2 = 17, respectively. These values are obtained for ISN2.0. For GSN, these numbers correspond to SN 1 = 130 and SN 2 = 10. By using the factor of 8.3 found for the emergence frequency of ARs between activity maximum and minimum by Harvey (1993) as a constraint, we obtain ∆m = 0.0946 (and thus m 1 = −2.5954 and m 2 = −2.7846).
Note that the values of m reach values higher than m 1 or lower than m 2 at times when the sunspot number is higher than SN 1 or lower than SN 2 , respectively. In particular, when the sunspot number is zero, the corresponding m values are m(0) = −3.952 for ISN2.0 and m(0) = −3.677 for GSN. The value of m(0) is different for ISN2.0 and GSN because the value of α, the free parameter of the model (see Table 1 Figure 2a shows the daily (black) and the monthly-smoothed (red) ISN2.0, while the evolution of m SN computed from the monthly-smoothed ISN for the solar cycle 21 is shown in Fig. 2b.  the analysis, whereas the five free parameters that are allowed to vary when fitting the data sets described in Sect. 2 are given in the lower part. One positive feature of the new model is that it has less free parameters than the old model (the old model required two free parameters to describe the emergence of ERs; see Table 1 of Wu et al. 2018a), giving it less 'wriggle room' for reproducing observational data. All parameters relevant to the emergence rates of ARs and SSEs (Eqs. 9-12) have been described in the previous section. Here we additionally comment on the decay and transfer times of the various components of the magnetic flux used in the ordinary differential equations describing the flux evolution (Eqs. 1-6).

Parameters of the model
The decay times of the ARs and SSEs, τ 0 AR and τ 0 SSE , are estimated using the observations by Parnell et al. (2009) and Thornton & Parnell (2011). Whereas Thornton & Parnell (2011) have analysed the emergence rate of different features as a function of their flux, Parnell et al. (2009) analysed the magnetic flux for all (i.e. instantaneously) observed features. By dividing the total number of the features of a given flux observed at a given instance by their emergence rate, we arrive at their mean lifetime. The lifetime of the features increases with their sizes or fluxes. For our purpose we make a simplification and calculate the lifetimes of ARs and SSEs as averages over all regions with fluxes above and below φ AR , respectively, thus obtaining τ 0 AR ≈ 10 days and τ 0 SSE ≈ 6 minutes. Since the features with the smallest flux dominate the power law distribution, τ 0 AR is short compared with lifetimes of large ARs and is closer to lifetimes of small ARs (e.g., Table 1 of the review by van Driel-Gesztelyi & Green 2015). For the same reason the SSE lifetime is close to that of internetwork elements, rather than of ERs. To compare better with observations of ERs, we therefore introduce φ ER = 10 18 Mx (see Fig. 1), which denotes roughly the lowest magnetic flux contained within ERs. If we now consider as ERs only the regions with φ ER < φ < φ AR , then we obtain a lifetime τ 0 ER ≈ 2 hours, which is comparable to the lifetimes of ≈hours to a day for regions with fluxes between 3 × 10 18 Mx and 10 20 Mx listed by van Driel-Gesztelyi & Green (2015). The maximum to minimum change of the flux emerging in ERs (i.e. with fluxes φ ER < φ < φ AR ) in our model is roughly a factor of 2.5, which is consistent with the results by Harvey (1993) and Harvey & Zwaan (1993). We stress that this distinction into ERs and internetwork fields is only used for comparison purposes. Within the model they are not distinguished (see Sect. 3.1).
The To do this, we use the genetic algorithm PIKAIA (Charbonneau 1995), which searches for the set of parameters minimising the difference between the modelled and the reference data sets. We minimise the sum of the reduced chi-squared values, χ 2 , taking the errors of the observations and the number of data points into account. In other words, we search for the maximum of 1/(χ 2 total + χ 2 OMF ); see  for details and Dasi-Espuig et al. (2016) for a discussion of uncertainties in the parameter fitting. The best-fit values of the parameters are listed in Table 1.
For comparison, the last column of Table 1 also lists the values of the five free parameters that are already present in the VS2010 model, as obtained by Wu et al. (2018a) for the most recent version of the model. This version used the same ISN2.0 input record extended back with the data from Vaquero et al. (2015) as done here. The values of the parameters are very close in the two models, except τ s SSE , for which Wu et al. (2018a) obtained 20.6 years compared to our 10.15 years (10.11 years when using GSN). Note, however, that Wu et al. (2018a) considered ERs rather than SSEs. Interestingly, the values we obtain are close to the value of 10.08 years found to produce a best-fit by  and is within the range 10-90 years that is consistent with observation (see . The value of about 4 years for the decay time of the slow open flux τ s open is close to the radial decay term with a timescale of 5 years introduced into surface flux transport simulations by Schrijver et al. (2002) and Baumann et al. (2006) to act on the unipolar fields at the solar poles. This was needed to reproduce the observed polar field reversals. Similarly, the τ r open obtained here (Table 1) is consistent with the estimate of the decay of closed flux carried by interplanetary coronal mass ejections, between 1 and 5 AU of between 38 and 55 days by Owens & Crooker (2006).

Reconstruction of the total and open magnetic flux
The computed total magnetic flux from 1965 onward is plotted in Fig. 3. Following Krivova et al. (2007), to account for the flux undetected due to the limited spatial resolution of observations (see Krivova & Solanki 2004) (Reimer et al. 2013). Note that the underlying 14 C data and thus also the OF recon- structed from them are decadal averages. The agreement between our model and the 14 C-based reconstruction is quite good. Particularly impressive is the agreement in the level of the open flux during the Maunder minimum, which is where we expect to see the biggest improvement relative to the old model. We emphasise that this 14 C-based record was not used to constrain our model. (Note that in the old model, the GSN record was used without the data from Vaquero et al. (2015) over the Maunder minimum, and the computed open flux in the old model is therefore essentially flat at the zero level.) Interestingly, over the 19th and the first half of the 18th centuries, the GSN-based reconstruction is closer to the isotope-based open flux than the reconstruction from the ISN, which lies somewhat higher.
A quantitative comparison of the total and open magnetic fluxes resulting from the old and the new models with the observations and independent reconstructions is presented in Table 2. The table lists the relative difference in means and the χ 2 values (in brackets) between the models and the data. For the total magnetic flux, the results are quite similar for both versions of the model. In both cases, the mean modelled total magnetic flux is somewhat closer to the observations when the GSN is used as input. The absolute difference in the means is slightly higher or lower for the new model if ISN2.0 or GSN are used, respectively. The new χ 2 values are somewhat higher than in the old model. This is, however, primarily due to fact that, by model design, the variability of the ER component on time scales shorter than the solar cycle in the old model was essentially smoothed out (see Eqs. 7-8) resulting in weaker short-term fluctuations than in the new model (see Fig. 4). Thus, if we smooth the total flux from both models with a 3-months window before comparing, the χ 2 values for the old model remain almost unchanged (0.036 for both ISN2.0 and GSN), while those for the new model decrease to 0.043 for ISN2.0 and 0.047 for GSN. In all cases, the χ 2 values are quite low.
For the open flux, the new model provides a notably better fit than the old model to all three alternative datasets. In all but one case, the absolute mean differences are significantly lower for the new model. The only exception is the GSN-based reconstruction versus the in-situ data by Owens et al. (2017), for which the absolute mean difference is slightly lower for the old model. However, the results are quite close for both versions of the model in this case. Note also that these data cover only a short recent period of time, over which the two models do not differ significantly. The χ 2 values are lower for the new model in all cases. The fit is poorest for the reconstruction based on the decadal values of 14 C. Very recently, new 14 C-based activity measures with annual resolution were published by Brehm et al. (2021). An application of our model to the dataset of Brehm et al. (2021) will be subject of a separate publication (Usoskin et al., submitted).  Lockwood et al. 2014), while the bottom part lists independent data sets that were not used for the optimisation (OF reconstruction from 14 C data by Wu et al. 2018b and the in-situ measurements by Owens et al. 2017). (*) F or decadally-averaged reconstructions. As 14 C data used for the reconstruction by Wu et al. (2018b) are decadal averages, only decadallyaveraged values of the OF could be reconstructed. Thus, to compute the corresponding χ 2 values, our reconstructions were re-sampled, too.

Summary and discussion
We have revised the simple model describing the evolution of the Sun's global total and open magnetic flux, originally proposed by Solanki et al. (2000Solanki et al. ( , 2002 and improved by . The new version of the model takes into account the observation that fluxes of magnetic features follow a single power law, including internetwork fields, ERs and ARs (Parnell et al. 2009;cf. Anusha et al. 2017). It also takes into account the fact that emergence rates of magnetic bipoles with fluxes between 10 16 Mx and 10 23 Mx, i.e. from the smallest ERs (and large internetwork features) to the largest ARs, also follow a power law according to the analysis by Thornton & Parnell (2011). We assume that the difference in emergence rates between the maximum and the minimum of solar activity is adequately described by the varying power-law exponent, affecting magnetic features with fluxes φ > 10 16 Mx. Thus, for the smallest features there is no change in emergence rate over the solar cycle, while the ratio of emergence rates during maximum to minimum increases steadily with increasing magnetic flux. Thus, the number of emerging ARs varies between the solar maximum and minimum of cycle 21 by a factor of 8.3, while this factor is 2.5 for ERs and close to unity for internetwork fields. These values are consistent with the respective ratios found by Harvey & Zwaan (1993) and Harvey (1993), while the fact that the flux in internetwork fields hardly changes over the cycle is in agreement with the results obtained by, e.g., Buehler et al. (2013) and Lites et al. (2014).
Using the sunspot number time series as input, the model returns time series of the open and total magnetic flux of the Sun. The main novel feature of the results of the model is that, in contrast to the earlier versions of the model, the output open magnetic flux does not drop to essentially zero during the Maunder minimum when almost no sunspots were present for multiple decades, in agreement with open flux reconstructed from 14 C data (e.g., Wu et al. 2018b). This significant improvement is a result of the model allowing for a non-zero emergence of magnetic flux in small-scale bipolar features (encompassing ERs and internetwork fields) even during extended periods when sunspots were not present, e.g. during the Maunder minimum and other, similar grand minima (e.g., Usoskin et al. 2007).
Even with this major update, the model still has some room for further improvements. Because it uses sunspot numbers as input (the only data of solar activity available prior to the middle of the 19th century), it cannot properly treat variations in solar activity that are not reflected in the number of sunspots. This is particularly evident during grand minima. During such times sunspots are only occasionally visible, whereas cosmogenic isotopes continue to display cyclic variations. This suggests that, in the context of the present model, the slope m continues to vary in a cyclical manner and can go lower than the lowest value we have obtained (m(0) = −3.952 for ISN2.0 and m(0) = −3.677 for GSN).
In the old version of the model the ER emergence was constructed as a smooth, sinusoidal function. In the new model, ER (and internetwork) emergence rate closely follows that of sunspots and has therefore the same temporal resolution as the input sunspot number series. As a consequence, the new model does not feature temporal lags or shifts to the corresponding sunspot cycle, as in the old model. However, it does account for the finding by Harvey (1994) that small-scale features (in her study ERs) belonging to a given cycle start emerging at a relatively high rate well before the sunspot cycle starts. The internetwork is independent of sunspot emergence and is, thus, not associated with a particular sunspot cycle. These regions are simply a result of the dynamo not having completely switched off even at times when there are no sunspots visible. In this way, the overlap between neighbouring solar magnetic cycles is naturally introduced. This was the main feature of the original model of Solanki et al. (2000) responsible for the change in the level of the Sun's open magnetic flux from one solar minimum to another, first noticed by Lockwood et al. (1999). At the same time, the new model allows for the dynamo to continue working and produce activity cycles during grand minima, which are sufficiently strong to modulate cosmic rays (see, e.g., Fig. 4 for the evolution of the open flux during the Maunder minimum) and hence influence the production of cosmogenic isotopes (Beer et al. 1998;Fligge et al. 1999;Usoskin et al. 2001), but are too weak to produce more than occasional sunspots.
The higher level of the magnetic flux during periods of very low solar activity (e.g. during the Maunder and the Dalton minima, see Fig. 4) will presumably lead to a weaker secular variation of the total solar irradiance (TSI) and in particular the rise of the TSI since the Maunder minimum (cf., Yeo et al. 2020), an important parameter for understanding the influence of solar irradiance variability on long-term climate trends (e.g., Gray et al. 2010;Solanki et al. 2013). The influence of the revised magnetic flux time series on TSI will be the subject of a future investigation.
Another important application of this model will be in reconstructing total magnetic flux and sunspot numbers from production rates of cosmogenic isotopes. Such an application will A&A proofs: manuscript no. AA_2021_40504 require the model to be inverted, as has successfully been done with the older version of the model (e.g., Usoskin et al. 2003Usoskin et al. , 2004Solanki et al. 2004;Usoskin et al. 2016;Wu et al. 2018b). Such work is ongoing and will be the subject of a follow-on publication.