Issue 
A&A
Volume 607, November 2017



Article Number  A76  
Number of page(s)  15  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201730689  
Published online  17 November 2017 
Parameter optimization for surface flux transport models
^{1} Department of Mathematical Sciences, Durham University, Durham DH1 3LE, UK
email: tim.j.whitbread@durham.ac.uk
^{2} Southwest Research Institute, 1050 Walnut St #300, Boulder, CO 80302, USA
email: amunozj@boulder.swri.edu
^{3} National Solar Observatory, Boulder, CO 80303, USA
email: gpetrie@nso.edu
Received: 23 February 2017
Accepted: 28 July 2017
Accurate prediction of solar activity calls for precise calibration of solar cycle models. Consequently we aim to find optimal parameters for models which describe the physical processes on the solar surface, which in turn act as proxies for what occurs in the interior and provide source terms for coronal models. We use a genetic algorithm to optimize surface flux transport models using National Solar Observatory (NSO) magnetogram data for Solar Cycle 23. This is applied to both a 1D model that inserts new magnetic flux in the form of idealized bipolar magnetic regions, and also to a 2D model that assimilates specific shapes of real active regions. The genetic algorithm searches for parameter sets (meridional flow speed and profile, supergranular diffusivity, initial magnetic field, and radial decay time) that produce the best fit between observed and simulated butterfly diagrams, weighted by a latitudedependent error structure which reflects uncertainty in observations. Due to the easily adaptable nature of the 2D model, the optimization process is repeated for Cycles 21, 22, and 24 in order to analyse cycletocycle variation of the optimal solution. We find that the ranges and optimal solutions for the various regimes are in reasonable agreement with results from the literature, both theoretical and observational. The optimal meridional flow profiles for each regime are almost entirely within observational bounds determined by magnetic feature tracking, with the 2D model being able to accommodate the mean observed profile more successfully. Differences between models appear to be important in deciding values for the diffusive and decay terms. In like fashion, differences in the behaviours of different solar cycles lead to contrasts in parameters defining the meridional flow and initial field strength.
Key words: magnetohydrodynamics (MHD) / Sun: activity / Sun: magnetic fields / Sun: photosphere
© ESO, 2017
1. Introduction
Surface flux transport (SFT) is a crucial component of the 11year sunspot cycle. SFT models have been developed and used for decades (e.g., DeVore et al. 1984; Wang et al. 1989a; van Ballegooijen et al. 1998; Schrijver & Title 2001; Baumann et al. 2004; Jiang et al. 2010) with some success, though results can be sensitive to the choice of parameters.
Surface flux transport models track the evolution of magnetic regions, which appear on the solar surface due to the rise of buoyant flux tubes (Fan 2009). Generally they emerge as bipolar magnetic regions (BMRs) with a leading polarity and a trailing polarity with respect to the eastwest direction. The leading polarities during a single cycle are usually the same across a hemisphere, with the reverse pattern occurring on the opposite hemisphere, following Hale’s polarity law (Hale 1924). Furthermore, the polarities reverse at the end of each cycle, giving rise to a magnetic cycle with a period of approximately 22 years. At the start of a cycle, magnetic regions emerge at latitudes of approximately ± 30°. As the cycle progresses, emergence regions migrate equatorwards before reaching ± 5° latitude at the end of the cycle. This emergence pattern was noted by Spörer (1880), and can be observed in the wellknown “butterfly diagram”.
Bipolar magnetic regions tend to emerge tilted at an angle with respect to the eastwest line (using the line connecting the centres of the opposing polarities), with the leading polarity emerging closer to the equator. This effect is due to the helical convective motions in the convection zone and is more pronounced at higher latitudes, according to Joy’s law (Howard 1991). In addition, the granular and supergranular convection cells provide a means of diffusion whereby flux is transported to the edge of the convection cells, spreading out across the solar surface and resulting in the leading flux cancelling across the equator with the corresponding opposite flux (Leighton 1964). The remaining trailing flux is transported to the polar regions via a combination of turbulent diffusion and a meridional flow (Howard 1979). Meridional circulation is a relatively slow transport mechanism, with speeds of ~10–20 m s^{1} observed via helioseismology by e.g., Braun & Fan (1998), Zhao & Kosovichev (2004), Jackiewicz et al. (2015). However, helioseismic recordings are near the limit of credibility (Komm et al. 2013) and so the flow profile is not well constrained. Flow speeds within this range have also been found by e.g., Komm et al. (1993) via the tracking of small magnetic features, and by Topka et al. (1982) via the comparison of polar zone filament distribution and polar magnetic field evolution. When the remaining flux eventually reaches the pole it cancels with the polar field from the previous cycle, culminating in polar field reversal. The reversal typically occurs at cycle maximum. For a historical review of surface flux transport, see Sheeley (2005).
There is evidence to suggest that the strength of the polar field at the end of the cycle is a good indicator of the strength of the following cycle (e.g., Schatten et al. 1978; MuñozJaramillo et al. 2013). Moreover, the strength of a cycle is anticorrelated with the duration of the cycle, according to the Waldmeier effect (Wolf & Brunner 1935). An explanation is offered by Cameron & Schüssler (2016) whereby activity belts near the base of the convection zone cancel across the equator with the opposing polarity. Stronger cycles have wider belts and so cancellation occurs earlier, resulting in a shorter cycle since all declining phases are approximately the same.
There is an everpresent need to predict the amplitude and length of future solar cycles since solar weather can have adverse effects on satellites, astronauts and technological systems on Earth. Given that disturbances such as coronal mass ejections and solar flares are usually attributed to regions of strong magnetic field, tracking and modelling magnetic regions on the solar surface, i.e., surface flux transport, has been highlighted as a key method for predicting and understanding solar cycle variability (e.g., Upton & Hathaway 2014b; Hathaway & Upton 2016), without using the more complicated calculations involved in dynamo simulations (for a review of dynamo theory, see Charbonneau 2014).
As mentioned above, the output of current SFT models is highly dependent on the parameters used. Parametrizations of the transport processes have been made, particularly for diffusion and meridional flow, but the exact forms are still debated, and are not necessarily in line with the limited observations available. Parameter studies of varying scope have been undertaken (e.g., Schrijver et al. 2002; Durrant & Wilson 2003; Baumann et al. 2004; Yeates 2014), but without complete parameter coverage. In this paper we aim to systematically find optimal parameters to be used in SFT models which accurately reproduce such features of the solar cycle as poleward flux transport “surges”, polar field reversal time, and polar field strength. The results can also be used to constrain the surface components of dynamo simulations to produce the most accurate cycle predictions to date. A similar study was performed by Lemerle et al. (2015), who used the same genetic algorithm used in this paper to find optimal parameters for a 2D SFT model for Cycle 21 only, with the view of coupling it to a 2D flux transport dynamo model (Lemerle & Charbonneau 2017). In contrast we analyse two distinct models with different dimensionality, namely 1D and 2D, and different dataassimilation techniques, initially for Cycle 23. While Lemerle et al. (2015) used the comprehensive BMR database compiled by Wang & Sheeley (1989), we use the same BMR database as Yeates et al. (2007) for the 1D model, and extract individual active regions from synoptic magnetograms for the 2D model. We also apply the 2D model to other cycles to search for cyclical variation.
In Sect. 2, we present the 1D model and the genetic algorithm used to perform the optimization, including a prescribed error structure dependent on latitude and magnetic field strength to factor in observational uncertainty. We also discuss the results of the optimization for Cycle 23. In Sect. 3, we present the 2D model which directly assimilates active regions into the simulation, and run the optimization process on this model for Cycle 23. In Sect. 4, we compare our optimal meridional flow profiles from both models with observations. Finally, we perform optimizations on the 2D model for Cycles 21, 22, and 24 in Sect. 5, before concluding in Sect. 6.
2. Onedimensional model
Since we aim to optimize the parameters of global, axisymmetric flows using the longitudeaveraged butterfly diagram, we start by considering a 1D SFT model. Using idealized BMRs as input, this model has the advantage of allowing many realizations to be run rapidly. In Sect. 3, we will compare the parameters selected by the 1D model to a 2D model with more realistic assimilation of active region data.
2.1. Numerical method
In the SFT model, the radial component of the magnetic field after averaging in longitude, B^{(}θ,t^{)}, evolves according to the advection–diffusion equation (Leighton 1964): (1)where R_{⊙} is solar radius, η is diffusivity, representing the diffusive effect of granular convective motions, τ is an exponential decay term added by Schrijver et al. (2002) to improve regular polar field reversal, and S is a source term for newly emerging magnetic regions. The profile describes poleward meridional flow. The magnetic field is decomposed into spherical harmonic form: (2)and the resulting coupled ordinary differential equations for the coefficients are solved using an adaptive RungeKutta (4, 5) timestepping method (Dormand & Prince 1980). The equations are solved on a grid of 180 cells equally spaced in latitude.
In the 1D model, new magnetic regions are assumed to take the form of idealized BMRs. Each BMR has a specified day of emergence; longitude and latitude; size; magnetic flux, including polarity; and tilt angle taken from an existing observational dataset where these properties were determined individually for each BMR from NSO synoptic magnetograms (Yeates et al. 2007). Using these data, the 1644 recorded BMRs from Cycle 23 (1st June 1996–3rd August 2008) are averaged in longitude and inserted into the model on the corresponding days of emergence. The BMR data are freely available at the Solar Dynamo Dataverse^{1} (Yeates 2016).
Because the meridional flow is a relatively slow transport mechanism, there is still a fair amount of uncertainty regarding its properties. A flexible profile is chosen to factor observational uncertainty into the optimization: (3)where increasing p produces a steeper gradient at low latitudes and a peak closer to the equator (Fig. 1). The amplitude v_{0} is chosen to be the maximum of  v , so increasing v_{0} increases the height of the peak. Both v_{0} and p are left as free parameters to be optimized. Lemerle et al. (2015) used a similar, but more sophisticated profile, which is discussed in Sect. 5.1. This provides substantially more flexibility, but introduces extra parameters into the optimization runs which could hinder convergence to a global maximum. A combination of exponential and sinusoidal functions adapted to helioseismic observations (Gizon & Duvall 2004) was utilized by Schüssler & Baumann (2006) for the meridional flow profile. Although this may better represent the actual meridional circulation, it cannot be accurately reproduced by the functional form in Eq. (3). In any case, the true functional form of the observed meridional circulation is uncertain, particularly at high latitudes.
Fig. 1 Three examples of the meridional flow profile in Eq. (3) for v_{0} = 15 m s^{1}, p = 2 (cyan), p = 5 (magenta) and p = 10 (black). 
In order to define the initial conditions, we use the profile of Svalgaard et al. (1978): (4)where B_{0} is left as a parameter to be optimized. Figure 2 (red) shows a crude fit of this expression to the observed profile from 1910 CR (blue). Whilst the observed profile is asymmetric across the equator and reveals some activity present at the equator, the prescribed profile represents a typical cycle minimum and ensures that the choice of initial profile is not hindering the optimization process, but rather aiding it with some flexibility. With this form of initial condition we can also be consistent across all regimes. Optimization runs performed using the observed initial condition show that the choice of initial condition in fact has a negligible effect on the optimal butterfly diagram.
2.2. Groundtruth data
As groundtruth data for optimization of the model, we use radialcomponent magnetogram data from US National Solar Observatory, Kitt Peak, in the form of fulldisk images. Prior to 2007 CR, these came from the Kitt Peak Vacuum Telescope, while from 2007 CR onwards we use Synoptic Optical Longterm Investigations of the Sun (SOLIS) data^{2}. To minimize noise in the polar regions of the map, we correct the butterfly diagram by calculating a cubic spline interpolation at each latitude of annual average measurements of highlatitude fields (poleward of ± 75°) which were observed with a preferable solar rotation axis tilt angle. A combination of real and interpolated data is used for the regions between ± 60° and ± 75° (see Petrie 2012). The resulting butterfly diagram is interpolated onto a uniform time grid at daily intervals. This is averaged over periods of 27 days, smoothed using a Gaussian filter to bring the unsigned flux down to a comparable level to the simulation, and finally sampled at the resolution of once per Carrington rotation. It should be noted that recently the SOLIS magnetograms were uniformly reprocessed and recalibrated. The butterfly diagram used here was made after this reprocessing.
2.3. Optimization algorithm
To search for optimal parameter sets where the model matches the observed butterfly diagram, we use the genetic algorithm PIKAIA 1.2. It was written by Charbonneau & Knapp (1995) at the High Altitude Observatory (HAO) of the National Center for Atmospheric Research (NCAR) and is publicly accessible^{3}.
PIKAIA is an evolutionary algorithm originally written in FORTRAN77. It is particularly effective for multimodal optimization problems, taking advantage of crossover and mutation operators which can induce jumps in parameter space, allowing for greater exploration and reducing the chances of getting trapped at a local maximum which is not the global maximum.
Fig. 2 Comparison between initial magnetogram (blue) and the profile given in Eq. (4) (red) with B_{0} = 8 G. 
The algorithm generates multiple random parameter sets, each entry within a userdefined range, and performs a model simulation for every parameter set, or “population member”. The population are ranked by “fitness” according to a userdefined “fitness function” which, in the case of model optimization, would usually be a comparison between a real reference case map and a modelgenerated map.
The highly ranked population members have a greater chance of being selected for the crossover or “breeding” process whereby sections of corresponding parameter strings from two members are interchanged to produce “offspring”, with the aim being that an individual will be produced with desirable features of both “parents” and become the fittest in the population. To increase variability and hence the likelihood for population improvement, random mutation is included, though this is a much slower process than crossover.
The crossovermutation process is run over a userdefined number of generations. Whilst PIKAIA is inherently stochastic and so finding an “acceptable” fit is never guaranteed, a large enough choice for the evolution period should ensure that the combined effect of the crossover and mutation operators has enough time to discover sufficiently fit population members.
For a more detailed introduction to PIKAIA and genetic algorithms in general, see Charbonneau (2002a,b).
In order to reduce the duration of computationally intensive optimizations, Metcalfe & Charbonneau (2003) created MPIKAIA, a freely accessible implementation of PIKAIA 1.2 in MPI^{4}. Rather than using a single processor for all model evaluations, a network of computers is used, and each of the model evaluations from a single generation is sent to a separate processor and computed simultaneously, achieving nearperfect parallelization. The time taken to complete the optimization therefore is entirely dependent on the number of generations. For example, for the 1D optimization problem in this paper, each model simulation takes approximately 90 s, so, with 46 processors, the runtime is brought down from 24 days to about 12.5 h for 46 population members evolved over 500 generations. These choices of population size and evolution period should be large enough to obtain good fits to the data, and will be used for the remainder of the paper unless specified.
2.4. Accounting for uncertainties in observations
The χ^{2} statistic is used as a measure of fit between the real and simulated butterfly diagrams: (5)where n is the number of gridpoints and x is the vector of k free parameters. Since improving best fit is a minimization process and PIKAIA is set up to maximize functions, the reciprocal of the measure, χ^{2}, is used as the necessary fitness function. The variance σ^{2} describes the error in both the measurements and the models, and we assume the form: (6)The variance plays two roles in the optimization. Firstly, it gives a meaningful value to the χ^{2} statistic. This allows us to compare the performance of different parameter combinations and time periods. Secondly, it effectively weights distinct locations (θ_{i},t_{j}) differently in the optimization, since the observed errors are assumed to have the form: (7)where ϵ is some small increment to ensure that the error is nonzero even in regions of low B_{obs}. This error structure reflects the uncertainties and inconsistencies in photospheric magnetic field observations (e.g., Riley et al. 2014). The factor of sinθ allows for the fact that the errors are in the original lineofsight measurements, whereas B_{obs} is the inferred radial field (Svalgaard et al. 1978). Overall, the effect of this error structure reduces the weight of observations both near the pole and in strong active regions. The resulting optimization will favour accuracy in the midlatitude “transport regions”. A further improvement of the χ^{2} statistic would be to include correlation between data points, but for simplicity here we assume independence.
Since the model error structure is unknown, we compute χ^{2} with σ_{model} = 0. This is sufficient for the purpose of comparing different model runs against the same set of observational data, with a higher value of χ^{2} indicating models that give a better match. The simulations are not sufficiently detailed to achieve a significant match at, say, the 99% level, which is evident from visual inspection of the butterfly diagrams. To achieve such a close match would be very challenging, since the large numbers of degrees of freedom n−k ~ 16 000–30 000 mean that the 99% interval for the χ^{2} statistic is narrow, typically [0.98,1.02].
In principle, we could estimate σ_{model} by increasing it and broadening the 99% interval until the value of χ^{2} falls within this interval. This would give a meaningful estimate of the “model error” in a particular run. But this would not change the ordering of different model runs, or indeed the final optimal parameters, so we have not included such analysis here.
2.5. Results
2.5.1. Optimal parameters
In addition to v_{0}, p and B_{0} as mentioned above, η and τ are also included in the optimization, resulting in 5 parameters to be optimized initially. Maximum and minimum limits are prescribed based on results from literature and observations (e.g., Schrijver et al. 2002; Hathaway & Rightmire 2010; Yeates 2014; Lemerle et al. 2015):

1.
100 km^{2} s^{1} ≤ η ≤ 1500 km^{2} s^{1};

2.
5 m s^{1} ≤ v_{0} ≤ 30 m s^{1};

3.
0 ≤ p ≤ 16;

4.
0 yr ≤ τ ≤ 32 yr;

5.
0 G ≤ B_{0} ≤ 50 G.
It should be noted that these ranges are deliberately made wider than results from literature to allow for a deeper exploration into the parameter space and to provide a better understanding of the SFT model. Table 1a shows the results of the optimization, with the corresponding optimal butterfly diagram in the top panel of Fig. 3, and the interpolated NSO data for Cycle 23 discussed in Sect. 2.2 in the bottom panel for direct qualitative comparison.
Optimal parameter sets for each optimization regime.
The equatorward migration of active regions is well represented by the BMR data, and large poleward surges are reproduced by the model. While the southern polar field reversal is well approximated by the model, the reversal in the northern hemisphere has a delay of approximately 5 CR. Furthermore, there are multiple weak poleward surges in the simulated butterfly diagram, most noticeably around 2000–2020 CR, which do not appear in the real butterfly diagram. This is likely to be a byproduct of approximating regions as BMRs and overestimating the contribution of flux from smaller regions. This buildup of flux results in a strong polar field that extends to lower latitudes, requiring a short decay timescale as is found in the optimization.
2.5.2. Parameter analysis
Fig. 3 Top: butterfly diagram for the optimal parameter 5set for the 1D model in Table 1a. Bottom: ground truth data for Cycle 23. 
During an optimization run, every single population member generated by PIKAIA can be recorded, and so a range of “acceptable” values can be obtained for each parameter. These can be found in square brackets below each optimal value in Table 1. The upper and lower bounds are taken to be the largest and smallest values for each parameter which produce fits above 95% of the maximum χ^{2}. Anything within these limits is classed as “acceptable”, though it must be noted that choosing to fix one parameter can alter the optimal solutions and bounds for others. Figure 4 shows such bounds, denoted by the left and right vertical purple lines on each plot, for all parameter populations from the optimization run that produced the optimal set in Table 1a. The optimal values are highlighted by the central vertical purple lines. Using the limits for v_{0} and p, acceptable meridional flow profiles were also found which are represented by the purple shading in the bottom right panel. The bold purple profile represents the optimum profile.
Fig. 4 Scatter plots of every population member for each parameter. The horizontal purple line denotes 95% of the maximum χ^{2}. The central vertical purple line is the optimum value for each parameter, with error bars given by the neighbouring vertical purple lines. The vertical blue lines in the bottom left and middle panels are the values obtained from fitting the velocity profile in Eq. (3) to observational data (see Sect. 4). The bottom right panel shows the optimal meridional flow profile (bold purple) with acceptable profiles represented by the surrounding purple shading. 
The diffusion parameter η has not yet been accurately measured, though some indirect measurements by Mosher (1977) and Komm et al. (1995) have found values within the range of 100–300 km^{2} s^{1}. Early simulations by Leighton (1964) used values up to 1000 km^{2} s^{1}, though studies by Baumann et al. (2004), Wang et al. (1989b) and Wang & Sheeley (1991) decreased it to ~ 600 km^{2} s^{1}, before Wang et al. (2002b) reduced it further to 500 km^{2} s^{1}. Our optimal value of 351.6 km^{2} s^{1}, however, is in better agreement with Yeates (2014), who found that η ∈ [200,450] km^{2} s^{1} produced a reasonable correlation between the butterfly diagrams, and Lemerle et al. (2015) who found an optimal value of 350 km^{2} s^{1} within an acceptable range of 240–660 km^{2} s^{1} for Cycle 21. Furthermore, Schrijver (2001) and Thibault et al. (2014) found diffusion coefficients of 300 km^{2} s^{1} and 416 km^{2} s^{1} respectively for randomwalkbased models, and Cameron et al. (2016) recently used a diffusion of 250 km^{2} s^{1}. The acceptable range in Table 1a is broad but can be attributed to multiple degrees of freedom in the optimization. The range covers most values discussed above.
The largescale meridional flow is poorly constrained by observations, as discussed in Sect. 2.1. Nevertheless, our optimal value of v_{0} = 14 m s^{1} is in accordance with both the observations and simulations. Doppler measurements by Ulrich (2010) estimated the maximum velocity to be between 14–16 m s^{1} for Cycles 22 and 23. Hathaway & Rightmire (2010) obtained an average maximum velocity of 10–12 m s^{1} for Cycle 23 via magnetic feature tracking, though crucially they observed that the flow is slower (approximately 8 m s^{1}) at cycle maximum and faster (11.5–13 m s^{1}) at minimum. This timedependence could be added to the model for greater realism, though it is not immediately clear how it could be parametrized in the optimization. Furthermore, they note that many SFT models use meridional flows which go to zero poleward of ± 75° latitude which is not necessarily what is observed, as well as other deviations from observations. Upton & Hathaway (2014a) prescribed a profile with a maximum velocity of 12 m s^{1} and Baumann et al. (2004) used 11 m s^{1}. Yeates (2014) discovered that a range of 11–15 m s^{1} improves butterfly diagram correlation, and Wang et al. (1989b) and Wang & Sheeley (1991) found that a range of 7–13 m s^{1} was acceptable. Wang et al. (2002b) found that a maximum velocity of 20–25 m s^{1} accurately reproduced solar cycle features, although they used a profile which differs significantly from observations.
The theta component of the flexible meridional flow profile in Eq. (3) was also used by MuñozJaramillo et al. (2009). They obtained a value of p = 2 by taking an average of helioseismic data weighted by density and fitting it to the sinusoidal profile. In this case p = 2 does not quite fall into the narrow acceptable range for p. The bottom right panel of Fig. 4 shows that values within this range generally correspond to a peak velocity at ± 30° before slowing down to 0 m s^{1} at ± 75°. As discussed above, this is not necessarily in line with observations. Taking every member of the population above the threshold, we find that the Pearson’s correlation coefficient between the acceptable values for v_{0} and p is r = 0.86, indicating that increasing the maximum velocity of the meridional profile generally requires an increase in p. A faster velocity means that active regions are transported away from the equator quicker. To counteract this, a larger value of p narrows the band of latitudes at which the velocity is fast, and additionally brings the maximum velocity closer to the equator.
Another interesting result is that of 2.4 yr for the exponential decay time τ. Schrijver et al. (2002) found that a decay time of 5–10 yr was necessary to replicate regular polar field reversal, and Yeates (2014) found that a decay time of 10 yr produced a better fit between real and simulated butterfly diagrams. Lemerle et al. (2015) found that exponential decay did not have a large effect on the polar field reversal and decided to set τ = 32 yr, effectively removing the decay term from the model. However, our optimal value for τ is close to the lower prescribed limit. This could be because of the model trying to account for the unusually weak polar field at the end of Cycle 23 while, for example, Lemerle et al. (2015) performed the optimization for Cycle 21.
Figure 5 highlights the need for the decay term in the 1D model when modelling Cycle 23. The purple curve, which represents the axial dipole moment obtained using the optimal parameter set in Table 1a, provides a better fit to the observed axial dipole moment (blue) than the peach curve, which is produced from the same parameter set but with the decay term omitted. In this case the polar field becomes too strong and is not weakened enough without the additional decay term. Jiang et al. (2015) found that the decay term was not required to obtain a close match between observed and simulated axial dipole moments, when using active region data from Li & Ulrich (2012). As well as using different active region data, a reduction in tilt angles and a smaller value of η = 250 km^{2} s^{1} were included to account for the lack of the decay term. If we use similar parameters for the 1D model, a better axial dipole moment fit is obtained at the expense of an accurate butterfly diagram. Hence we stress that the optimal values in Table 1 are with respect to the measure of choice in Eq. (5), and other choices of metric might give different results.
Of course, the choice of decay term is not independent of the other parameters, and the Pearson’s correlation coefficient between the acceptable values of v_{0} and τ is r = 0.81: an increase in the flow speed corresponds to less trailing flux being transported to the poles, so a fast decay to weaken the polar field would not be required in the presence of a faster flow.
It should be noted that the decay term in Eq. (1) is not directly observed and was added by Schrijver et al. (2002) to produce regular polar field reversals over multiple cycles. Wang et al. (2002a) overcame this problem by increasing the meridional flow speed for stronger cycles. Baumann et al. (2006) gave a physical explanation of the decay term; namely, it is the effect of radial (i.e., inward) diffusion of flux into the solar interior, which cannot be accounted for directly in the SFT model. In spherical harmonics, different modes decay at different rates, whereas in the exponential decay used by Schrijver et al. (2002), all modes decay at the same rate. Baumann et al. (2006) found that the lowestorder mode decayed the slowest at a rate of 5 yr (with a corresponding volume diffusion of η = 100 km^{2} s^{1}), in good agreement with the findings of Schrijver et al. (2002). When we include this more sophisticated form of radial diffusion in our model and perform the optimization, we find the lowestorder mode to have an optimal decay time of τ_{1} = 2.7 yr (with a corresponding volume diffusion of η = 190 km^{2} s^{1}), in good agreement with the decay time found in Table 1a. Because of this good agreement, we opt to continue to use the original exponential decay parameter.
The optimal value for B_{0} is significantly higher than that used to approximate the initial profile in Fig. 2. This might be attributed to the choice of functional form in Eq. (4); not enough flux is prescribed between ± 45° and ± 80°, so the algorithm compensates for this by increasing the maximum flux at ± 90°. Alternatively, a strong initial polar field is also required to counteract the short decay time needed to reproduce the weak polar field at the end of Cycle 23.
Fig. 5 Axial dipole moments calculated from observed data (blue), the parameter set in Table 1a (purple), and the same parameter set but with the decay term omitted (peach). 
2.6. Tilt angles
Fig. 6 Top: butterfly diagram for the optimal parameter 6set for the 1D model with reduced tilt angles in Table 1b. Bottom: ground truth data for Cycle 23. 
Some studies (e.g., Yeates 2014; Jiang et al. 2011) found that multiplying the tilt angle of each BMR by a scaling factor reduces the polar field strength and improves polar field reversal, since the reduced tilt inhibits equatorial crosscancellation and hence each magnetic region will contribute less to the axial dipole. A multiplicative tilt angle factor (TAF) is included here as an extra parameter to be optimized within the range 0 ≤ TAF ≤ 1.5. Table 1b shows the results for the 6 parameter case, with the corresponding butterfly diagram in the top panel of Fig. 6.
The optimal value of 0.55 for TAF is lower than that found by Yeates (2014) (TAF ~ 0.8) and Jiang et al. (2011) (TAF ~ 0.72). It predictably produces a weaker polar field than in the case above where it was not included. Given that the main aim of the algorithm is to reduce differences between the real and simulated butterfly diagram pixels, it is reasonable to expect that the optimization algorithm will rely heavily on diffusion and high amplitudes of meridional flow to achieve weak polar fields, although it should be noted that this effect is reduced by the weighting in σ. Introducing the tilt angle factor as a means of reducing the polar field allows for the decay time to increase and the maximum meridional flow velocity to decrease, suggesting a delicate balance between the parameters and the roles they play in the model. While the polar field strength is better approximated in this case, the active regions are much weaker than in the 5parameter case, and polar field reversal occurs later in the simulation in both hemispheres. The fitness value of χ^{2} = 1.09 is above the 99% interval given in Sect. 2.3, which seems to indicate that the model matches the observations better than a randomly chosen map from the observed distribution. This is plainly a limitation of the χ^{2} statistic; in particular, it likely indicates the presence of a significant σ_{model} term possessing a more complex structure over θ and t. In principle, it could be caused by too large a prescribed σ_{obs}, or by the relatively strong assumption of independence, or possibly overfitting of the model, although the latter is unlikely given the small number of parameters in the model. It should be noted that the scaling of tilt angles is not a physical phenomenon, rather a method of reducing the flux in the model, though Cameron et al. (2010) argued that scaling the tilt angle by a factor of 0.7 mimics the effect of inflows around active regions. Moreover, DasiEspuig et al. (2010) found an inverse correlation between cycle strength and tilt angle, suggesting that tilt angle variation plays a significant role in polar field variation.
3. Twodimensional model
Yeates et al. (2015) developed a 2D model^{5} which assimilates specific shapes of magnetic regions into the simulation on the day of emergence. The aim of the model is to better assimilate strong, multipolar regions, which are not accurately portrayed in a simpler bipolar form, as in the 1D model above, with the hope of simulating a more realistic photospheric field. This selection feature requires the model to be 2D. The model is fully automated, providing consistent highlighting of strong magnetic regions, and is designed to replace preexisting regions rather than superimposing new ones. The SFT equation for the radial component of the magnetic field in 2D, B(θ,φ,t), is: (8)where ω represents differential rotation and all other parameters are identical to those in Eq. (1). Note that differential rotation averaged out and so played no role in the axisymmetric 1D model. Rather than using a spectral method like the 1D model, the evolution equations (for the vector potential) are solved in the Carrington frame using a finitedifference method on a spatial grid of 180 cells equally spaced in sinelatitude and 180 cells equally spaced in longitude. Unlike meridional flow, differential rotation is well constrained by observations, and in the model is parametrized as (Snodgrass & Ulrich 1990): (9) The 2D model contains a parameter BPAR which determines the threshold above which magnetic flux is assimilated into the simulation, in the form of individual strongflux regions. Yeates et al. (2015) chose the threshold of BPAR = 15 G in order that the difference between the observed unsigned flux and simulated unsigned flux (due to the smoother magnetic field distribution) remained approximately constant. This parameter is subsequently added to the optimization. If given enough freedom, the algorithm would gradually reduce BPAR, allowing more and more magnetic regions to be inserted until the original synoptic map is essentially copied in (analogous to BPAR ~ 0 G). To avoid this, the lower bound is set at 10 G with an upper bound of 50 G. Figure 7 shows snapshots of 1928 CR from four simulations with alternative values of BPAR between 10 G and 50 G, and all other parameters fixed. As the threshold BPAR increases, fewer active regions are assimilated into the simulation.
Fig. 7 Four snapshots of 1928 CR from simulations with active regions selected by different values of the magnetic flux threshold BPAR and all other parameters fixed. Here BPAR increases from left to right and top to bottom. 
Fig. 8 Top: butterfly diagram for the optimal parameter 5set for the 2D model with varying BPAR in Table 1c. Bottom: ground truth data for Cycle 23. 
3.1. Fiveparameter optimization
The synoptic magnetograms from NSO Kitt Peak are used to identify strong regions for assimilation. For simplicity, Yeates et al. (2015) did not incorporate exponential decay into the model as in Eq. (1). We perform optimization runs for the model both without decay and with the decay term included. Initially we consider the former case. Aside from BPAR, parameters are given the same upper and lower limits as in Sect. 2.5.1. Table 1c shows the results of the optimization. The corresponding butterfly diagram is shown in the top panel of Fig. 8.
The 2D model qualitatively improves the butterfly diagram, with active regions predictably more accurate, leading to the inclusion of more poleward surges in the simulation which can be identified in the observed butterfly diagram (though the gradient and strength of each surge is not always correct), and a more realistic polar field. The optimal parameters in Table 1c are within the range of other results from simulations and observations described in Sect. 2.5.2. A diffusivity of η = 455.6 km^{2} s^{1} is a stronger diffusivity than in the 1D model, but the inclusion of an exponential decay term is expected to reduce this. An increased diffusivity is somewhat supported by Virtanen et al. (2017), who used a value of η = 400 km^{2} s^{1} in the same 2D model but for a single simulation of multiple cycles. The range and optimal value for v_{0} is lower than for the original 1D case, indicating that there can be inherent differences between models. Moreover, Virtanen et al. (2017) found that a value of v_{0} = 11 m s^{1} correctly reproduced shapes of poleward surges and polar fields, in excellent agreement with our optimal value.
Figure 9 shows every generated value of BPAR against χ^{2}. The central vertical line indicates the optimum value of 39.8 G, with the left and right vertical lines denoting the acceptable range for BPAR, as in Fig. 4. The value of 15 G used by Yeates et al. (2015) is outside of this range, and for the remainder of the 2D optimizations, BPAR is fixed at the optimal value of 39.8 G to attain consistency. This should ensure that only newly emerging regions are inserted for each Carrington rotation. However, the presence of the strong midlatitudinal region of positive flux in the northern hemisphere at 2000–2020 CR could be attributed to the choice of large BPAR, since smaller regions of negative flux which would otherwise cancel out this positive flux are not being assimilated. The bottom left panel of Fig. 7 closely represents the scenario when BPAR is set at its optimal value. Virtanen et al. (2017) used a threshold of BPAR = 50 G, and this lies just outside of our acceptable range. Comparing the bottom two panels of Fig. 7, however, shows that the differences between our optimal value and their chosen value are minor.
Fig. 9 Each population member for the 5parameter optimization of the 2D model with BPAR restricted to BPAR ≥ 10 G. The horizontal green line denotes 95% of the maximum χ^{2}. The central vertical green line is the optimum value for each parameter, with error bars given by the neighbouring vertical green lines. 
3.2. Incorporating exponential decay
As discussed above, the decay parameter τ was originally added to the SFT model to produce regular polar field reversals. The 2D model did not initially take account of this decay time, but we incorporate it to assess whether the optimal value in Table 1a is reasonable.
As shown in Fig. 10, including the decay term improves timing of polar field reversal by 5–10 CR, but is not enough to replicate the observed reversal time. Poleward surges are generally wider in the simulation, leading to the reduction of some midlatitude features, most notably the strong surge of positive flux at 2000–2020 CR in the northern hemisphere, which is more visible in Fig. 8.
The optimization results are shown in Table 1d. Surprisingly, the addition of an extra decay term induces a minimal decline in diffusion, and it is not enough to bring it down to 351.6 km^{2} s^{1} as found in the 1D case. Rather, B_{0} increases to account for the stronger decay of the polar fields in this regime. Most significantly, we obtain an optimal value of τ = 4.5 yr. This is higher than the optimum found in the 1D model and in closer agreement with Schrijver et al. (2002), although the acceptable range is considerably wider towards the upper limit, indicating that a decay term may not be required in the assimilative model. This is supported by the value of χ^{2} which does not increase significantly with the addition of the decay term. Furthermore, Fig. 11 shows that the axial dipole moments calculated using the optimal parameter sets for the 2D model, with and without the exponential decay term (brown and green curves respectively), both produce good fits to the observed profile (blue). This indicates that the method of new flux assimilation in the 2D model is better able to account for the weak polar field at the Cycle 23/24 minimum than the idealized BMRs used in the 1D model, since it does not require an additional decay term. Coupled to the short optimal decay timescale are smaller optimal values for v_{0} and p, suggesting that the relationships and correlations discussed in Sect. 2.5 also hold for the 2D case.
Fig. 10 Top: butterfly diagram for the optimal parameter 5set for the 2D model with varying τ in Table 1d. Bottom: ground truth data for Cycle 23. 
Fig. 11 Axial dipole moments calculated from observed data (blue), the parameter set in Table 1c (green), and the parameter set in Table 1d (brown). 
4. Comparison with meridional flow observations
Although observations of the meridional flow are not yet fully reliable, we can use the data that are available to try to add a further constraint to the optimization.
David Hathaway kindly provided us with measurements of the meridional flow for Solar Cycle 23, calculated by tracking features in images from the Solar and Heliospheric Observatory (SOHO), Michelson Doppler Imager (MDI; Scherrer et al. 1995). The data were supplied as coefficients of the following parametrization: (10)The meridional flow measurements for each Carrington rotation are shown in Fig. 12 (blue curves). The observations tend to follow either a fast or slow flow, highlighted by denser blue areas, indicating the dependence on time and that the flow transitions between the two extremes throughout the cycle. Additionally, for a small number of Carrington rotations an equatorward counterflow is observed at high latitudes, though it should be noted that such a counterflow was not visible in HMI data (Hathaway & Upton 2014). The choice of flexible profile in Eq. (3) does not allow for this phenomenon.
The optimal profile using the parameters from the 1D optimization in Table 1a is shown in purple in Fig. 12 for comparison. Whilst the observed and optimal profiles are similar in shape, the optimal profile is too fast and reaches its peak at a slightly lower latitude. Moreover, the observed profiles tend to extend beyond ± 75° but the optimal profile chooses to go to zero throughout the polar regions, giving a possible explanation as to why many SFT models incorporate this feature. Furthermore, the 1D optimal profile remains almost completely within the bounds given by the observations, excluding at its peak in the northern hemisphere for which asymmetry in the observations can be held responsible.
The green and brown profiles in Fig. 12 represent the optima for the 2D model excluding and including exponential decay respectively. Both profiles are fully contained within the observational limits, except for a small section of the brown curve in the southern hemisphere which is due to a lower than average maximum velocity. Of the three optimal profiles, the 2D regime without decay matches the average observed profile the closest, whilst the decayenhanced flow is slightly slower (though Hathaway & Rightmire 2010, observed speeds of 8 m s^{1} at cycle maximum). It does, however, continue to latitudes poleward of ± 70°, almost emulating the observational data. One limitation of tracking magnetic features to measure the meridional flow is that it is not always easy to distinguish between the effects of the meridional flow and the effects of supergranular diffusion. For this reason, flows derived from feature tracking tend to peak at higher latitudes (e.g., Dikpati et al. 2010, Fig. 1), giving a possible explanation as to why the observed curves in Fig. 12 tend to peak at higher latitudes than the modelled curves.
Fig. 12 Comparison of various meridional flow profiles: observed for each CR (blue), 1D optimum (purple), 2D optimum (green) and 2D optimum with decay (brown). 
We use a nonlinear leastsquares fitting method to fit the parametrized form of the meridional flow in Eq. (3) to the average observed coefficients given by David Hathaway to ensure it is actually possible to match the observed profile. The average observed and fitted profiles, shown in Fig. 13 (blue and red respectively), match closely for v_{0} = 11.3 m s^{1} and p = 1.87, and slight asymmetry in the average observed profile is confirmed. This value of p is close to that of MuñozJaramillo et al. (2009) and is within the acceptable ranges for p in the above 2D regimes, but is outside the equivalent range in the 1D optimization run, whence we infer that the 1D model requires the maximum velocity to be closer to the equator than is observed.
Fig. 13 Comparison of average observed (blue) and fitted (red) meridional flow profiles. 
Given that the parametrization is able to closely fit the observed data, we could fix one of the velocityrelated parameters, say p, to the observed value and perform optimization runs for the two models. We choose p because the model is generally less sensitive to the choice of v_{0}, and p = 1.87 is outside the acceptable range for the 1D model.
The optimization results with p fixed in the 1D model are shown in Table 1e. The value p = 1.87 corresponds to a maximum velocity at ± 35°, meaning poleward transport is slower at low latitudes. This results in more flux cancellation across the equator and so more trailing flux is present in the transport regions, as observed in the top panel of Fig. 14. This feature appears to be a common occurrence in the standard SFT model (cf. Figs. 3 and 6). The upshot of this numerically is that the selected decay time of 1.9 yr is even shorter than in the original 1D case to counteract the large amounts of flux accumulating at the poles. This couples with a slow velocity, made even slower by the small value of p, adhering to the relationship found in Sect. 2.5. The timing of polar field reversal, meanwhile, is reproduced reasonably accurately. Except for a marginally smaller value of χ^{2}, fixing p does not significantly hinder the quantitative performance of the 1D model, even though p = 1.87 is not in the acceptable parameter range for regime (a).
Fig. 14 Top: butterfly diagram for the optimal parameter 4set for the 1D model with fixed p = 1.87 in Table 1e. Bottom: ground truth data for Cycle 23. 
With the higherlatitudinal velocity peak and the absence of τ in the 2D model, the resulting diffusion value given in Table 1f is slightly larger than in previous regimes. Contrary to expectation, the optimal maximum velocity is higher than the previous 2D cases, but still with wide error bounds. Given that p = 1.87 lies within the acceptable range in regime (c), it is reasonable to expect that optimal values and associated ranges would be in line with results in Sect. 3 and hence observations and previous studies. Consequently the optimal butterfly diagram (top panel of Fig. 15) confirms this, offering only subtle changes to Fig. 8, for example a polar field restricted to higher latitudes due to the increase in diffusivity.
Fig. 15 Top: butterfly diagram for the optimal parameter 3set for the 2D model with fixed p = 1.87 in Table 1f. Bottom: ground truth data for Cycle 23. 
5. Other solar cycles
With its automated assimilation of active region data, the 2D model can easily be adapted for other cycles, provided there is sufficient data available. Evaluations of Cycles 21, 22 and 24 (up to the end of 2015) using NSO data have been carried out to search for cycletocycle variation.
5.1. Cycle 21
Table 1g shows the optimum parameters for Cycle 21 (1st May 1976−10th March 1986). Both η and v_{0} are in agreement with previous studies. Most notably, v_{0} = 9.2 m s^{1} is slower than the maximum speed of Cycle 23, supporting Upton & Hathaway (2014a): a faster flow in Cycle 23 would have resulted in a weaker polar field at cycle minimum since leading flux is taken away from the equator quickly and so has less time to cancel across the equator. This optimum value, however, is just outside the range of 10–13.2 m s^{1} as found by Komm et al. (1993) using feature tracking during Cycle 21. Conversely, this range overlaps with a large portion of the 95% confidence interval obtained by the optimization population.
The interpolated NSO data is shown in the bottom panel of Fig. 16 with the corresponding simulated butterfly diagram in the top panel of Fig. 16. Aside from a negativepolarity observational artefact in the northern hemisphere at 1680 CR, many features of active regions are well reproduced. There are three instances of large concentrations of opposite flux being transported polewards in the northern hemisphere; the latter of these is overestimated by the simulation and this could be attributed to the model incorrectly reading in the corresponding emergence region. Polar field reversal for both poles is too late in the model, particularly in the northern hemisphere where the difference is in the region of 10 CR.
Fig. 16 Top: butterfly diagram for the optimal parameter 4set for the 2D model in Table 1g. Bottom: ground truth data for Cycle 21. 
Lemerle et al. (2015) performed a similar optimization process for Cycle 21 using a 2D model and a BMR database compiled by Wang & Sheeley (1989). Although they used a different parametrization for the meridional flow and different sources of flux, their optimal parameter ranges for η and v_{0} were in good agreement with those in Table 1g. The diffusion coefficient η = 455.7 km^{2} s^{1} lies within their acceptable range of 240–660 km^{2} s^{1} and v_{0} = 9.2 m s^{1} falls between 8–18 m s^{1} as calculated by PIKAIA in their study. They used the following functional form to represent meridional flow: (11)The optimization returned values of v_{0} = 12 m s^{1}, q = 7, v = 2 and w = 8. This gave a profile similar to that of Wang et al. (2002b), but with a less extreme steep gradient at the equator. However, when normalized, the profile shape was comparable to the observed profile formed from Doppler measurements obtained by Ulrich (2010), and the observed profile lay well within the error bars for the optimal solution, except for some return flows at high latitudes, which were not incorporable in Eq. (11), mirroring the limitation of our parametrization in Eq. (3). Using a nonlinear leastsquares fitting method, we are able to attempt to fit the functional form in Eq. (3) to the versatile meridional profile in Eq. (11). The best fit corresponds to values of v_{0} = 13.6 m s^{1} and p = 3.88. This value for v_{0} is in agreement with observations and acceptable ranges for other regimes, but is above the range for Cycle 21. Despite lying within the acceptable range, p = 3.88 favours the high values for p obtained from optimization runs as opposed to the lower values extracted from observational data. This could suggest an inherent flaw within the SFT model whereby the model performs better when the maximum velocity is prescribed to be closer to the equator.
5.2. Cycle 22
Table 1h shows the optimization results for Cycle 22 (10th March 1986–1st June 1996). The fit is marginally worse than for Cycle 21, but optimal values for η and v_{0} remain within in plausible ranges. The optimal diffusion in this case increases to 506.2 km^{2} s^{1}, but is in better agreement with Wang et al. (2002b). The optimal maximum velocity for Cycle 22 is even smaller than that of Cycle 21, further supporting the fact that a slower meridional flow results in a stronger polar field at cycle minimum, and explaining the high optimal maximum velocity for Cycle 23. van Ballegooijen et al. (1998) performed SFT simulations for Cycle 22 with η = 450 km^{2} s^{1} and v_{0} = 11 m s^{1} which produced polar field strength in agreement with observations. Again, these values are in accordance with ranges given in Table 1h.
The ground truth data is shown in the bottom panel of Fig. 17 and the simulated butterfly diagram is in the top panel of Fig. 17. The model has recreated polarity reversal much more successfully here, with only a slight delay in the north. Towards the end of the cycle there is a large buildup of positive flux and some weak, but visible, poleward surges in the northern hemisphere that have appeared in the simulation but are not observed in the real butterfly diagram.
Fig. 17 Top: butterfly diagram for the optimal parameter 4set for the 2D model in Table 1h. Bottom: ground truth data for Cycle 22. 
5.3. Cycle 24 (so far...)
Table 1i shows the results for the first half of Cycle 24 (3rd August 2008–1st January 2016). We obtain a much higher value of χ^{2} for Cycle 24 compared to previous cycles, but we suspect that this is due to the relative ease of modelling only half a cycle as opposed to modelling longterm effects. The diffusivity η = 454.6 km^{2} s^{1} is within viable ranges found in literature, though the maximum velocity is close to the lower prescribed bound. The initial polar field B_{0} = 4.2 G is lower than in previous cycles as the model needs to replicate the weak polar field at the Cycle 23/24 minimum. Acceptable ranges of parameters are generally broad, but performing the optimization on the full cycle in the next few years should tighten the upper and lower bounds. Indeed, when a similar optimization process is performed on half of Cycle 23, the acceptable ranges are found to be wider, though the shorter time period has a negligible effect on the specific optimal values.
The interpolated Kitt Peak data is shown in the bottom panel of Fig. 18 with the corresponding simulated butterfly diagram in the top panel of Fig. 18. Although a large portion of the cycle is yet to take place, there are still some notable features, such as the prominent leadingpolarity region between 2100 CR and 2110 CR in the northern hemisphere. This region was the primary subject of Yeates et al. (2015). Polar field reversal is slightly late in the simulated butterfly diagram; performing an optimization once the full cycle has completed might remedy this, though a region of negative polarity in the northern hemisphere at 2160 CR may not correctly be reproduced, unless the data is corrected.
Fig. 18 Top: butterfly diagram for the optimal parameter 4set for the 2D model in Table 1i. Bottom: ground truth data for Cycle 24. 
Including exponential decay in the model for Cycles 21, 22 and 24 produces optimal values of τ = 10.2 ∈ [3.1,32.0] yr, τ = 7.6 ∈ [3.1,32.0] yr and τ = 15.1 ∈ [2.5,32.0] yr respectively. These are in better agreement with Schrijver et al. (2002) and Lemerle et al. (2015), indicating that the low optimal value for τ may only be necessary in Cycle 23 in order to successfully reconstruct the unusually weak polar fields at Cycle 23/24 minimum.
6. Conclusions
The aim of this paper was to use a genetic algorithm to find optimal parameters to be used in surface flux transport simulations, subsequently helping us understand the behaviour and interplay of the many physical processes on the Sun. We began by obtaining optimized parameter sets for a 1D SFT model for Cycle 23, both with and without a multiplicative tilt angle factor. From these simulations we obtained viable ranges for parameters. We found that these ranges and optimal solutions were in good agreement with results from previous studies and from observations. We also looked at the interaction of parameters, highlighting the positive correlations between the meridional velocity parameters v_{0} and p, and exponential decay time τ.
We repeated the optimization process on a 2D assimilative model and found that optimum parameters were mostly within ranges of those from the 1D case, but distinct enough to suggest that the differences between models could be important. We also found an optimum value for the assimilation threshold, which was significantly greater than used previously by Yeates et al. (2015). Qualitatively, the 2D model produced a more accurate butterfly diagram than the 1D model, particularly at the poles. We also included an exponential decay term in the 2D model which produced an optimal value of 4.5 yr, which lies outside the acceptable range found in the 1D case and is in agreement with the values obtained by other authors. Including decay induced a decrease in the velocity parameters, but given that the acceptable range extended to the upper limits of exploration, its inclusion may not be necessary in the 2D model. There is the possibility that we did not model decay realistically, which could have led to a strong polar field. That the 2D model was able to give an acceptable match to the observed butterfly diagram and axial dipole moment without a decay term is evidence that it is superior to the 1D model, which was unable to do so with the corresponding optimal parameters. It suggests that the method of flux assimilation in the 2D model is superior to the insertion of idealized BMRs, as used both in the 1D model and in most other SFT models.
We were then able to compare the optimal meridional profiles from different regimes with observations made from feature tracking. The profiles from regimes (a), (c), and (d) were each almost completely within the range of observed flows, but the 1D optimal profile was faster than the average observed flow, while the 2D profile with decay included was too slow. The 2D profile without an extra decay term, however, best matched the average observed profile. Fixing the observed profile in both models resulted in varied success; the 2D model was able to accommodate the observations comfortably, whilst the 1D model saw a reduction in most parameters and a butterfly diagram containing an excess of flux in the transport regions.
Finally, the optimization process was repeated using the 2D model for Cycles 21, 22, and 24, producing plausible results for Cycles 21 and 22; Cycle 24 may need more time to progress to capture the longterm effects of the cycle in the optimal parameters, particularly in narrowing some of the range of viable solutions, although an optimization run performed over the same period of time for Cycle 23 showed that the optimal parameters themselves are barely affected; it was just the ranges of acceptable values which widened due to fewer constraints. In order to predict the axial dipole moment at the Cycle 24/25 minimum and hence the amplitude and length of Cycle 25, randomly generated magnetic regions with properties based on empirical relations must be used to simulate the remainder of the cycle (e.g., Upton & Hathaway 2014b; Cameron et al. 2016). Analysis of multiple cycles highlighted significant differences in meridional circulation speed, supporting the evidence for slower meridional flows during stronger cycles, and initial profile strength, supporting the proposed relationship between cycle strength and polar field strength at the preceding cycle minimum. Our multiple cycle analysis also highlighted cycledependence of the decay term τ. At present, the best form and magnitude of such a decay term remain to be determined by the community. However, our results (and the others mentioned) do suggest that it can help to improve the match with observations, at least for Cycle 23. It is intriguing that it seems to be less important for the preceding cycles. This could either be because the decay is compensating for some other deficiency of the model that has changed in Cycle 23, such as the inability to reproduce the unusually weak polar field at the end of the cycle, or the radial diffusion of flux did really change from one cycle to the next, presumably due to some difference in the flows and magnetic field in the convection zone. This is an interesting subject for future study, but is beyond the scope of this paper where we consider only the surface. All optimization runs were performed with respect to a prescribed variance which was proportional to both latitude and B_{obs}, chosen to correspond to uncertainty in observational data. It should be noted that comparing fitness values is always with respect to the chosen error structure in this paper. For other studies of modelling Cycle 23 see, e.g., Schrijver & Liu (2008), Yeates et al. (2010), Yeates & Mackay (2012), Jiang et al. (2013).
While the flexibility in the problem is beneficial in the respect that it allows more freedom, it can also have drawbacks. For example, the choice of fitness function is crucial to deciding which regime or parameter choice is “best” for each model, but depends entirely on what the user regards as important. Lemerle et al. (2015) used a combination of χ^{2} statistics which measured the differences between real and simulated timelatitude maps, axial dipoles and “transport regions” (latitudes ± 34° to ± 54°). These statistics were balanced equally in the final fitness function. Weighting could have been applied in favour of particular features, though it is not obvious how best to put this into practice. Alternatively, weighting could be applied to different sections of the map, i.e., active, transport and polar regions, to force the algorithm to return parameters which produce those specific regions more accurately. We chose a comparison between the real and simulated timelatitude maps, with an associated error structure, as we considered the general reproduction of the whole map to be foremost in importance.
The adaptability of the 2D model provides a wide scope of possible future directions. One such direction is testing variability between different measuring instruments to ascertain whether inconsistent literature results could simply be due to the choice of observatory or satellite. This comes with the issue of either deciding on or computing an appropriate value for the assimilation threshold BPAR for different datasets. Another future possibility that takes advantage of the model’s assimilation technique is to optimize multiple cycles at the same time. We have shown that there exists variation in parameters between cycles, so a single optimal parameter set for more than one cycle would be unrealistic. An alternative method would be to treat each cycle separately, coupled only at each cycle minimum, where the final profile of the previous cycle becomes the initial profile of the next.
Our methodology assumes a static meridional flow. The inclusion of a timevarying meridional flow in the optimization could significantly alter results, however parametrizing timedependence without introducing too many parameters is not a trivial procedure. On the other hand, largescale inflows towards active regions were first observed by Gizon et al. (2001), and Cameron & Schüssler (2012) proposed that these flows were at least partially responsible for variation of meridional flow over the solar cycle. Indeed, MartinBelda & Cameron (2016) found that the inflows increased the effect of flux cancellation and also reduced the latitudinal separation of polarities, thereby decreasing the axial dipole moment contribution of a bipolar region. This process weakens the polar field in the same way that a timedependent meridional flow can, and although we have not accounted for inflows in this study, it is an option under consideration for future work. An alternative method for reducing the polar field is using a fluxdependent diffusion parameter whereby the presence of a strong magnetic field quenches diffusion (e.g., MuñozJaramillo et al. 2011).
In the near future we hope to use PIKAIA to optimize a kinematic 3D dynamo model (Yeates & MuñozJaramillo 2013) using the results in this paper to constrain the surface evolution.
Acknowledgments
T.W. thanks Durham University Department of Mathematical Sciences for funding his Ph.D. studentship. A.R.Y. thanks STFC for financial support. A.R.Y. and G.J.D.P. thank ISSI, Bern for supporting their collaboration as part of an International Team on global magnetic modelling. A.M.J. is very grateful to George Fisher and Stuart Bale for their support at the University of California, Berkeley, and Phil Scherrer for his support at Stanford University. This research is partly funded by the NASA Grand Challenge grant NNH13ZDA001N and NASA LWS grant NNX16AB77G. Data were acquired by SOLIS instruments operated by NISP/NSO/AURA/NSF. We are grateful to David Hathaway for also supplying data for this study and to Ian Vernon for his invaluable statistical advice as well as to the anonymous referee for useful suggestions.
References
 Baumann, I., Schmitt, D., Schüssler, M., & Solanki, S. K. 2004, A&A, 426, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baumann, I., Schmitt, D., & Schüssler, M. 2006, A&A, 446, 307 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Braun, D. C., & Fan, Y. 1998, ApJ, 508, L105 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R. H., & Schüssler, M. 2012, A&A, 548, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cameron, R. H., & Schüssler, M. 2016, A&A, 591, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cameron, R. H., Jiang, J., Schmitt, D., & Schüssler, M. 2010, ApJ, 719, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R. H., Jiang, J., & Schüssler, M. 2016, ApJ, 823, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P. 2002a, NCAR Tech. Note, NCAR/TN450+IA (Boulder, CO: National Center for Atmospheric Research), 1 [Google Scholar]
 Charbonneau, P. 2002b, NCAR Tech. Note, NCAR/TN451+STR (Boulder, CO: National Center for Atmospheric Research), 1 [Google Scholar]
 Charbonneau, P. 2014, ARA&A, 52, 251 [Google Scholar]
 Charbonneau, P., & Knapp, B. 1995, NCAR Tech. Note, NCAR/TN418+IA (Boulder, CO: National Center for Atmospheric Research), 1 [Google Scholar]
 DasiEspuig, M., Solanki, S. K., Krivova, N. A., Cameron, R., & Peñuela, T. 2010, A&A, 518, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DeVore, C. R., Boris, J. P., & Sheeley, Jr., N. R. 1984, Sol. Phys., 92, 1 [Google Scholar]
 Dikpati, M., Gilman, P. A., & Ulrich, R. K. 2010, ApJ, 722, 774 [CrossRef] [Google Scholar]
 Dormand, J. R., & Prince, P. J. 1980, J. Comput. Appl. Math., 6, 19 [CrossRef] [MathSciNet] [Google Scholar]
 Durrant, C. J., & Wilson, P. R. 2003, Sol. Phys., 214, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2009, Liv. Rev. Sol. Phys., 6, 4 [Google Scholar]
 Gizon, L., & Duvall, T. L. 2004, in MultiWavelength Investigations of Solar Activity, eds. A. V. Stepanov, E. E. Benevolenskaya, & A. G. Kosovichev, IAU Symp., 223, 41 [Google Scholar]
 Gizon, L., Duvall, Jr., T. L., & Larsen, R. M. 2001, in Recent Insights into the Physics of the Sun and Heliosphere: Highlights from SOHO and Other Space Missions, eds. P. Brekke, B. Fleck, & J. B. Gurman, IAU Symp., 203, 189 [Google Scholar]
 Hale, G. E. 1924, Proc. Nat. Acad. Sci., 10, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., & Rightmire, L. 2010, Science, 327, 1350 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., & Upton, L. 2014, J. Geophys. Res. (Space Phys.), 119, 3316 [Google Scholar]
 Hathaway, D. H., & Upton, L. A. 2016, J. Geophys. Res. (Space Phys.), 121, 10 [Google Scholar]
 Howard, R. 1979, ApJ, 228, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Howard, R. F. 1991, Sol. Phys., 136, 251 [Google Scholar]
 Jackiewicz, J., Serebryanskiy, A., & Kholikov, S. 2015, ApJ, 805, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, J., Işik, E., Cameron, R. H., Schmitt, D., & Schüssler, M. 2010, ApJ, 717, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, J., Cameron, R. H., Schmitt, D., & Schüssler, M. 2011, A&A, 528, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jiang, J., Cameron, R. H., Schmitt, D., & Schüssler, M. 2013, Space Sci. Rev., 176, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, J., Cameron, R. H., & Schüssler, M. 2015, ApJ, 808, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Komm, R. W., Howard, R. F., & Harvey, J. W. 1993, Sol. Phys., 147, 207 [Google Scholar]
 Komm, R. W., Howard, R. F., & Harvey, J. W. 1995, Sol. Phys., 158, 213 [NASA ADS] [Google Scholar]
 Komm, R., González Hernández, I., Hill, F., et al. 2013, Sol. Phys., 287, 85 [CrossRef] [Google Scholar]
 Leighton, R. B. 1964, ApJ, 140, 1547 [Google Scholar]
 Lemerle, A., & Charbonneau, P. 2017, ApJ, 834, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Lemerle, A., Charbonneau, P., & CarignanDugas, A. 2015, ApJ, 810, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Li, J., & Ulrich, R. K. 2012, ApJ, 758, 115 [NASA ADS] [CrossRef] [Google Scholar]
 MartinBelda, D., & Cameron, R. H. 2016, A&A, 586, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Metcalfe, T. S., & Charbonneau, P. 2003, J. Comput. Phys., 185, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Mosher, J. M. 1977, Ph.D. Thesis, California Institute of Technology, Pasadena [Google Scholar]
 MuñozJaramillo, A., Nandy, D., & Martens, P. C. H. 2009, ApJ, 698, 461 [NASA ADS] [CrossRef] [Google Scholar]
 MuñozJaramillo, A., Nandy, D., & Martens, P. C. H. 2011, ApJ, 727, L23 [NASA ADS] [CrossRef] [Google Scholar]
 MuñozJaramillo, A., DasiEspuig, M., Balmaceda, L. A., & DeLuca, E. E. 2013, ApJ, 767, L25 [Google Scholar]
 Petrie, G. J. D. 2012, Sol. Phys., 281, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Riley, P., BenNun, M., Linker, J. A., et al. 2014, Sol. Phys., 289, 769 [Google Scholar]
 Schatten, K. H., Scherrer, P. H., Svalgaard, L., & Wilcox, J. M. 1978, Geophys. Res. Lett., 5, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129 [Google Scholar]
 Schrijver, C. J. 2001, ApJ, 547, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., & Liu, Y. 2008, Sol. Phys., 252, 19 [Google Scholar]
 Schrijver, C. J., & Title, A. M. 2001, ApJ, 551, 1099 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., De Rosa, M. L., & Title, A. M. 2002, ApJ, 577, 1006 [NASA ADS] [CrossRef] [Google Scholar]
 Schüssler, M., & Baumann, I. 2006, A&A, 459, 945 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sheeley, Jr., N. R. 2005, Liv. Rev. Sol. Phys., 2, 5 [Google Scholar]
 Snodgrass, H. B., & Ulrich, R. K. 1990, ApJ, 351, 309 [Google Scholar]
 Spörer, G. 1880, Publikationen des Astrophysikalischen Observatoriums zu Potsdam, 2, 1 [Google Scholar]
 Svalgaard, L., Duvall, Jr., T. L., & Scherrer, P. H. 1978, Sol. Phys., 58, 225 [Google Scholar]
 Thibault, K., Charbonneau, P., & Béland, M. 2014, ApJ, 796, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Topka, K., Moore, R., Labonte, B. J., & Howard, R. 1982, Sol. Phys., 79, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Ulrich, R. K. 2010, ApJ, 725, 658 [NASA ADS] [CrossRef] [Google Scholar]
 Upton, L., & Hathaway, D. H. 2014a, ApJ, 792, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Upton, L., & Hathaway, D. H. 2014b, ApJ, 780, 5 [NASA ADS] [CrossRef] [Google Scholar]
 van Ballegooijen, A. A., Cartledge, N. P., & Priest, E. R. 1998, ApJ, 501, 866 [NASA ADS] [CrossRef] [Google Scholar]
 Virtanen, I. O. I., Virtanen, I. I., Pevtsov, A. A., Yeates, A. R., & Mursula, K. 2017, A&A, 604, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, Y.M., & Sheeley, Jr., N. R. 1989, Sol. Phys., 124, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y.M., & Sheeley, Jr., N. R. 1991, ApJ, 375, 761 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y.M., Nash, A. G., & Sheeley, Jr., N. R. 1989a, ApJ, 347, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y.M., Nash, A. G., & Sheeley, Jr., N. R. 1989b, Science, 245, 712 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wang, Y.M., Lean, J., & Sheeley, Jr., N. R. 2002a, ApJ, 577, L53 [Google Scholar]
 Wang, Y.M., Sheeley, Jr., N. R., & Lean, J. 2002b, ApJ, 580, 1188 [NASA ADS] [CrossRef] [Google Scholar]
 Wolf, R., & Brunner, W. 1935, Astronomische Mitteilungen der Eidgenössischen Sternwarte Zurich, 14, 105 [Google Scholar]
 Yeates, A. R. 2014, Sol. Phys., 289, 631 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R. 2016, Bipolar magnetic regions determined from NSO synoptic carrington maps, Harvard Dataverse [Google Scholar]
 Yeates, A. R., & Mackay, D. H. 2012, ApJ, 753, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., & MuñozJaramillo, A. 2013, MNRAS, 436, 3366 [Google Scholar]
 Yeates, A. R., Mackay, D. H., & van Ballegooijen, A. A. 2007, Sol. Phys., 245, 87 [Google Scholar]
 Yeates, A. R., Mackay, D. H., van Ballegooijen, A. A., & Constable, J. A. 2010, J. Geophys. Res. (Space Phys.), 115, A09112 [Google Scholar]
 Yeates, A. R., Baker, D., & van DrielGesztelyi, L. 2015, Sol. Phys., 290, 3189 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, J., & Kosovichev, A. G. 2004, ApJ, 603, 776 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Three examples of the meridional flow profile in Eq. (3) for v_{0} = 15 m s^{1}, p = 2 (cyan), p = 5 (magenta) and p = 10 (black). 

In the text 
Fig. 2 Comparison between initial magnetogram (blue) and the profile given in Eq. (4) (red) with B_{0} = 8 G. 

In the text 
Fig. 3 Top: butterfly diagram for the optimal parameter 5set for the 1D model in Table 1a. Bottom: ground truth data for Cycle 23. 

In the text 
Fig. 4 Scatter plots of every population member for each parameter. The horizontal purple line denotes 95% of the maximum χ^{2}. The central vertical purple line is the optimum value for each parameter, with error bars given by the neighbouring vertical purple lines. The vertical blue lines in the bottom left and middle panels are the values obtained from fitting the velocity profile in Eq. (3) to observational data (see Sect. 4). The bottom right panel shows the optimal meridional flow profile (bold purple) with acceptable profiles represented by the surrounding purple shading. 

In the text 
Fig. 5 Axial dipole moments calculated from observed data (blue), the parameter set in Table 1a (purple), and the same parameter set but with the decay term omitted (peach). 

In the text 
Fig. 6 Top: butterfly diagram for the optimal parameter 6set for the 1D model with reduced tilt angles in Table 1b. Bottom: ground truth data for Cycle 23. 

In the text 
Fig. 7 Four snapshots of 1928 CR from simulations with active regions selected by different values of the magnetic flux threshold BPAR and all other parameters fixed. Here BPAR increases from left to right and top to bottom. 

In the text 
Fig. 8 Top: butterfly diagram for the optimal parameter 5set for the 2D model with varying BPAR in Table 1c. Bottom: ground truth data for Cycle 23. 

In the text 
Fig. 9 Each population member for the 5parameter optimization of the 2D model with BPAR restricted to BPAR ≥ 10 G. The horizontal green line denotes 95% of the maximum χ^{2}. The central vertical green line is the optimum value for each parameter, with error bars given by the neighbouring vertical green lines. 

In the text 
Fig. 10 Top: butterfly diagram for the optimal parameter 5set for the 2D model with varying τ in Table 1d. Bottom: ground truth data for Cycle 23. 

In the text 
Fig. 11 Axial dipole moments calculated from observed data (blue), the parameter set in Table 1c (green), and the parameter set in Table 1d (brown). 

In the text 
Fig. 12 Comparison of various meridional flow profiles: observed for each CR (blue), 1D optimum (purple), 2D optimum (green) and 2D optimum with decay (brown). 

In the text 
Fig. 13 Comparison of average observed (blue) and fitted (red) meridional flow profiles. 

In the text 
Fig. 14 Top: butterfly diagram for the optimal parameter 4set for the 1D model with fixed p = 1.87 in Table 1e. Bottom: ground truth data for Cycle 23. 

In the text 
Fig. 15 Top: butterfly diagram for the optimal parameter 3set for the 2D model with fixed p = 1.87 in Table 1f. Bottom: ground truth data for Cycle 23. 

In the text 
Fig. 16 Top: butterfly diagram for the optimal parameter 4set for the 2D model in Table 1g. Bottom: ground truth data for Cycle 21. 

In the text 
Fig. 17 Top: butterfly diagram for the optimal parameter 4set for the 2D model in Table 1h. Bottom: ground truth data for Cycle 22. 

In the text 
Fig. 18 Top: butterfly diagram for the optimal parameter 4set for the 2D model in Table 1i. Bottom: ground truth data for Cycle 24. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.