© ESO 2019

1. Introduction

A major goal of solar physics is to understand the solar cycle: the near-periodic rise and decline of solar magnetic activity. Periods of maximum activity correspond to more frequent space weather events like solar flares and coronal mass ejections, which pose threats to satellites and astronauts and can cause technological disruption at Earth. On the other hand, the state of the Sun’s magnetic field at solar minimum gives us an indication of the general global behaviour of the next solar cycle (Schatten et al. 1978; Muñoz-Jaramillo et al. 2013; Pesnell 2016). It is generally accepted that solar magnetic activity is maintained by a dynamo mechanism. However, it is not yet possible to directly measure the magnetic fields in the solar interior. Instead we currently rely on mathematical models to provide insight into the dynamo process.

There are numerous varieties of dynamo model, each having their own strengths and limitations (for reviews see Charbonneau 2010, 2014). Here, however, we focus on Babcock–Leighton (B–L) models (Babcock 1961; Leighton 1964, 1969). In the B–L regime, toroidal field is converted to poloidal field via the emergence and decay of active regions at the photosphere. The cross-equatorial cancellation of leading polarity flux and subsequent preferential polewards transport of trailing flux results in a polarity reversal of the polar field. The polar field is then pumped down into the convection zone, where it is sheared back into toroidal field by differential rotation and transported equatorwards by the returning branch of meridional circulation or latitudinal turbulent magnetic pumping, becoming the seed field for the next cycle. An appealing property of B–L dynamos is that they operate in line with observations of the photospheric radial magnetic field (Wang et al. 1989; Wang & Sheeley 1991). Furthermore, they have been found to reproduce features of the solar cycle (Dikpati et al. 2004); Mackay & Yeates (2012).

Progression in this area thus far has primarily been through the implementation of 2D or 2  ×  2D B–L dynamo models (e.g. Wang et al. 1991; Durney 1995; Chatterjee et al. 2004; Guerrero & de Gouveia Dal Pino 2008; Lemerle & Charbonneau 2017; Bhowmik & Nandy 2018). However, we would ideally like to develop 3D B–L dynamo models in order to realistically model the emergence of buoyant magnetic structures and fully describe the evolution of magnetic fields under the effects of diffusion, differential rotation, and meridional circulation. These models are more complex and require in-depth calibration in order to match the observed magnetic field. Nevertheless, success in overcoming these obstacles would be a sizeable step towards the development of a forecasting model for the Sun-Earth system (Nita et al. 2018). This would hopefully provide us with the most accurate solar cycle predictions to date.

Yeates & Muñoz-Jaramillo (2013) developed KD3, a 3D kinematic B–L dynamo model. In KD3, the magnetohydrodynamic induction equation describes the evolution of the mean magnetic field:

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\mathit{v}\times \mathit{B})-\mathrm{\nabla }\times (\eta \mathrm{\nabla }\times \mathit{B}),\end{array}$$(1)

for a prescribed velocity field v(r,θ,ϕ,t) and prescribed turbulent diffusivity η(r). There is no small-scale α-effect. Equation (1) is solved in a spherical shell using a finite volume scheme. For more details see Appendix A of Yeates & Muñoz-Jaramillo (2013).

Unlike previous 2D B–L models, KD3 explicitly models the buoyant emergence of flux tubes through the convection zone (Fan 2009). In the 2D models, the active region emergence process has either been parametrised through a volumetric α-effect term in the induction equation, or through manual deposition of regions at the surface, corresponding to areas of strong toroidal field at the base of the convection zone (e.g. Durney 1997; Nandy & Choudhuri 2001; Muñoz-Jaramillo et al. 2010; Guerrero et al. 2012). The deposition method has also been used in another 3D B–L dynamo model, developed by Miesch & Dikpati (2014; see also Miesch & Teweldebirhan 2016; Karak & Miesch 2017, 2018; Hazra et al. 2017; Hazra & Miesch 2018). However, these “non-local” methods make magnetic flux conservation difficult to enforce because the process of forming the emerging region from the pre-existing toroidal field is not followed explicitly through the induction Eq. (1).

In KD3, a time-dependent velocity perturbation is included which is intended to capture the effects of advection and buoyancy on the flux tubes. A similar emergence method has been used by Kumar et al. (2018) and Kumar et al. (2019). The non-axisymmetric perturbation has a radial component, which transports the tube outwards through the convection zone to the surface; a vortical component, which models the helical convective motions and gives rise to tilts in the active regions; and a diverging component, responsible for expanding the tube as the density decreases. The tube centre velocity is set so that the travel time from r = 0.7 R_⊙ to r = R_⊙ is 25 days, after which the perturbation is removed. This method ensures the conservation of magnetic flux during the emergence process. This is in contrast to the deposition method, where active regions are disconnected and an interior structure is assumed only close to the surface.

Although the KD3 emergence approach is flux-conserving, and Yeates & Muñoz-Jaramillo (2013) showed that the model is able to reproduce the qualitative behaviour of active region decay at the surface, leading to polewards transport of flux and reversal of the polar field, closer inspection has shown that the quantitative details of the surface evolution are significantly different from 2D surface flux transport (SFT) models, even when the same horizontal flows and diffusivity and same initial B_r are used at the surface. As an example, the SFT evolution of a single bipolar magnetic region (BMR), shown in Fig. 1 and placed at 10° latitude with flux 1 × 10²² Mx and a tilt angle of 30°, is shown in the top panel of Fig. 2, and the KD3 equivalent is shown below. The parameters used are the same as those given in Sect. 3. The BMR is inserted in the SFT simulation at the time when the flux has stopped emerging in KD3, that is, when the unsigned flux at the photosphere has reached its peak (Fig. 3). Even though the differential rotation, meridional flow and horizontal diffusion in the SFT model match the surface parameters of the KD3 simulation, the transport to the poles is faster in the SFT case. In addition, the top panel of Fig. 3 shows that there is significantly more flux present at the surface in the KD3 system. There is also a large difference between the respective evolutions of the polar flux (bottom panel of Fig. 3). In KD3, the south polar field barely develops by the end of the simulation, and the peak of the north polar field is stronger and occurs three years later than in the SFT case.

Fig. 1. Three-dimensional image of an emerged active region in KD3. Magnetic field lines are connected to the toroidal field at the base of the convection zone and the radial magnetic field is shown at the transparent surface.

Fig. 2. Top: longitude-averaged evolution of Br for a single BMR in a 2D SFT model. Bottom: surface component of the 3D dynamo model showing the equivalent evolution of the same BMR.

Fig. 3. Top: comparison of unsigned surface flux from the 2D SFT simulation (blue) and 3D dynamo simulation (orange). Bottom: comparison of northern (solid and dotted lines) and southern (dashed and dash-dotted) polar flux from the same two simulations, where polar flux is defined as the flux polewards of 70° latitude.

Given that the SFT model parameters have been carefully calibrated to match the evolution of B_r on the real Sun (Lemerle et al. 2015; Whitbread et al. 2017), we start from the premise that it is the KD3 model that needs to be modified. In this paper, we show that the incorrect evolution of surface flux in the KD3 model arises from the fact that BMRs remain connected to the base of the convection zone for several years after emergence, owing to the low diffusivity in the convection zone. This effect is illustrated in Sect. 2 using a simplified Cartesian 2D model, where we show that increasing the convection zone diffusivity can improve the surface evolution. In Sect. 3, we verify that the same is true in KD3 when simulating a single region. However, Sect. 4 shows that increasing the diffusivity in a full-cycle simulation of KD3 has a catastrophic effect on the dynamo. The implications for future 3D modelling are discussed in Sect. 5.

2. 2D model of active region decay

We begin by investigating a 2D model that illustrates the basic cause of the difference between the KD3 and SFT models. Inspired by van Ballegooijen (1998), we take a 2D Ω-loop representing a newly-emerged BMR in the convection zone and evolve it according to diffusion alone. The benefit of a simpler toy model is that it captures the diffusive effects of a 3D model but is computationally less expensive, and at this stage we are not interested in other features such as the amount of poloidal field produced.

Here we use Cartesian co-ordinates (x,z) which denote the width and depth of the convection zone domain respectively, with −0.4 ≤ x ≤ 0.4 and 0.6 ≤ z ≤ 1. Neglecting variation in the y-direction, we write B in terms of a flux function as:

$$\begin{array}{c}\hfill \mathit{B}=\mathrm{\nabla }\times (A{\mathit{e}}_y).\end{array}$$(2)

Neglecting advection, Eq. (1) reduces to

$$\begin{array}{c}\hfill \frac{\partial A}{\partial t}=\eta \left(z\right){\mathrm{\nabla }}^{2}A.\end{array}$$(3)

Importantly, we allow the diffusivity η to be a function of z, so that we can investigate the effect of different diffusivity profiles with depth. The effect of advection will be considered in the 3D simulations of Sect. 3. We also simultaneously evolve a 1D surface diffusion model as the analogue of the SFT model. For visualization a potential-field extrapolation is performed in the corona. For our first initial condition, the region is assumed to have emerged and is connected to the toroidal field at the base of the convection zone (see left-hand panel of Fig. 4), as in KD3, and is of the form:

$$\begin{array}{c}\hfill A_0=exp(-\frac{z-0.6}{0.04})+\frac{1}{2}exp(\frac{{(z-1)}^2}{0.4}-\frac{x^2}{0.008}).\end{array}$$(4)

Fig. 4. Two initial conditions used in this paper: an active region connected to the toroidal field (left) and an active region disconnected from the toroidal field (right).

We impose periodic boundary conditions in x and set ∂A/∂t = 0 at the base (z = 0.6). At the surface (z = 1) we follow van Ballegooijen (1998) and van Ballegooijen & Mackay (2007) by setting:

$$\begin{array}{c}\hfill B_{x,\mathrm{c}z}=\beta B_{x,\mathrm{c}or},\end{array}$$(5)

where B_x,cz is the horizontal field at the convection zone boundary and B_x,cor is the horizontal component of a potential extrapolation into the corona. Then the parameter β determines whether the interior field at the photosphere is matched to the potential field in the corona (β = 1), or whether it is purely radial (β = 0), which was the original boundary condition in KD3 (and, indeed, in most other models). This will allow us to assess the effect of the top boundary condition on radial diffusion, although for most tests we set β = 0.

In general we will use the following depth-dependent two-step profile for η(z):

$$\begin{array}{cc}\hfill \eta \left(z\right)=& \frac{1}{{\eta }_{\max }}[{\eta }_{\mathrm{c}}+\frac{{\eta }_0-{\eta }_{\mathrm{c}}}{2}(1+\mathrm{erf}\left(\frac{z-R_1}{{\mathrm{\Delta }}_1}\right))\hfill \\ \hfill & +\frac{{\eta }_{\mathrm{s}}-{\eta }_0-{\eta }_{\mathrm{c}}}{2}(1+\mathrm{erf}\left(\frac{z-R_2}{{\mathrm{\Delta }}_2}\right))],\hfill \end{array}$$(6)

where η_max is chosen such that the maximum value of the diffusivity is 1. Here η_c is the core diffusivity, η₀ is the diffusivity in the convection zone, and η_s is the surface diffusivity. The step locations and thicknesses are R_i and Δ_i respectively. The profiles used in this paper are shown in Fig. 5. The solid orange curve shows the profile used by Yeates & Muñoz-Jaramillo (2013) in KD3, with parameters η_c = 10⁸ cm² s⁻¹, η₀ = 1.6 × 10¹¹ cm² s⁻¹, η_s = 6 × 10¹² cm² s⁻¹, R₁ = 0.71, Δ₁ = 0.03, R₂ = 0.95 and Δ₂ = 0.025. The other profiles will be described later in the section. For a given diffusion profile and boundary condition, Eq. (3) is solved using an explicit finite difference method with Euler timestepping.

Fig. 5. Normalised multi-step diffusion profiles used in this paper, against a log-scale. The solid orange curve is from KD3, the dashed purple curve is the profile which takes into account diffusivity quenching, and the dotted yellow curve is derived from mixing-length theory.

When the KD3 diffusion profile is used, it is clear from the top right panel of Fig. 6 that there is significantly more flux at the surface than would be expected without radial derivatives, as we saw for the KD3 model in Sect. 1. This is because the relatively low diffusion below z = 0.9 does not allow for much diffusive transport, and field lines remain attached to the toroidal field at the base of the convection zone. Because the field lines are fixed in place, movement at the surface is heavily restricted and cancellation at the boundary is limited, resulting in an excess of surface flux. This transpires even though diffusion is stronger near the surface, as indicated by the outwards bulging of field lines.

Fig. 6. Snapshots of magnetic field lines from simulations using the KD3 (top row), quenching (second row), MLT (third row) and constant (bottom row) diffusion profiles. For the first three profiles, the left-hand column shows the case where the region is connected to the base, and the middle column shows the simulation with a disconnected initial condition. The black dashed line is the top of the domain, above which is shown a potential-field extrapolation. Right-hand column: evolution of magnetic flux at the surface, compared to a 1D surface model (green dashed line). The bottom row (constant η) compares the effects of potential (blue) and radial (red) boundary conditions but using the same connected initial condition.

To demonstrate the different evolution for a disconnected active region such as considered by Miesch & Dikpati (2014; hereafter STABLE), we alter the initial condition slightly from Eq. (4):

$$\begin{array}{c}\hfill A_0=\frac{1}{2}exp(-\frac{x^2+{(z-1)}^2}{0.008}).\end{array}$$(7)

This forms a potential field below the surface, disconnected completely from the base of the convection zone (see right-hand panel of Fig. 4).

The top middle panel of Fig. 6 shows a snapshot from the simulation with the original KD3 diffusion profile and disconnected initial condition. Because field lines are no longer connected to the toroidal field at the base of the domain, the weak diffusion no longer plays a role in anchoring field lines in place. This allows for more diffusive transport and cancellation of magnetic flux at the surface. Instead of flux cancellation occurring at the side boundary after field lines are pushed outwards as before, cancellation takes place between the two polarities of the active region. The consequence is that there is less surface flux in the early stages of evolution in comparison to the 1D model. The cancellation rate eventually decreases, but we observe in the top right panel that there is less surface flux present at the end of the simulation than in the case where the connected initial condition was used. The upshot is that the disconnected region qualitatively provides a better match to the surface than the connected region.

In the presence of strong magnetic fields, turbulent diffusivity can be suppressed (Roberts & Soward 1975). This “quenching” can be included in models via a non-linear relationship whereby the diffusion parameter η is scaled by the reciprocal of the square of the magnetic field (e.g. Tobias 1996; Gilman & Rempel 2005; Muñoz-Jaramillo et al. 2008; Guerrero et al. 2009). By instead taking the geometric spatiotemporal average over many effective diffusivity profiles, Muñoz-Jaramillo et al. (2011) approximated the effect of the dynamically quenched diffusion using a fixed profile in the form of Eq. (6) by applying the following parameters: η_c = 10⁸ cm² s⁻¹, η₀ = 1.6 × 10¹¹ cm² s⁻¹, η_s = 3.25 × 10¹² cm² s⁻¹, R₁ = 0.71, Δ₁ = 0.017, R₂ = 0.895 and Δ₂ = 0.051. This is shown as the dashed purple curve in Fig. 5, and will henceforth be referred to as the “quenching profile” for simplicity.

A snapshot from the simulation using the quenching profile is displayed in the second row of Fig. 6. When the original connected initial condition is prescribed, the field lines diffuse downwards initially, but approximately halfway through the simulation the direction of motion changes and the magnetic field starts to diffuse upwards. We note a reduction in the surface flux, presumably because the stronger diffusivity levels extend deeper into the domain and the field lines have more freedom to move, allowing for more diffusive transport. However, we find again that flux cancellation is hindered by the weak diffusion in the lower convection zone, which keeps the field lines attached to the toroidal field.

If instead we start with a disconnected region (middle panel), we find that, as for the KD3 profile, flux cancels inwardly because field lines are not connected to the base of the convection zone. However, it diffuses at a much faster rate than the regime with the KD3 diffusion profile (and hence the 1D case), and by the end of the simulation the majority of the surface flux has been cancelled.

The third profile we experiment with is derived from mixing-length theory (MLT; Prandtl 1925). Muñoz-Jaramillo et al. (2011) used the solar interior model of Christensen-Dalsgaard et al. (1996) to estimate the mixing-length parameter α_p and hence the diffusivity profile based on GONG data. The value of diffusion found for the convection zone is up to two orders of magnitude larger than those used in KD3 and other kinematic dynamo simulations in literature. This is because simulated dynamo action has not yet been achieved in flux transport dynamos with such strong diffusion. Muñoz-Jaramillo et al. (2011) attempted to reconcile the MLT estimates with numerical values by incorporating diffusivity quenching, leading to the quenching profile above. Nevertheless, a fit to the MLT profile was also made in the form of Eq. (6), with the following resulting parameters: η_c  =  10⁸ cm² s⁻¹, η₀  =  1.4 × 10¹³ cm² s⁻¹, η_s  =  10¹⁰ cm² s⁻¹, R₁  =  0.71, Δ₁  =  0.015, R₂  =  0.96 and Δ₂  =  0.09. This profile is the dotted yellow curve shown in Fig. 5.

A snapshot from the corresponding simulation is shown in the third row of Fig. 6. With this diffusion profile and connected initial condition, the field initially diffuses downwards before being pushed back up due to the diffusion gradient at the surface. This surface flux then diffuses to the boundary where it cancels. Low diffusivity at the base means the field still remains attached to the toroidal field but a much larger diffusivity throughout the convection zone helps transport flux upwards from as deep as z = 0.7. We see that field lines are being pushed together at the top of the domain due to the reduced diffusivity near the surface and a balance between outwards and inwards diffusion. At a higher cadence, we observe that this causes the field lines to reconnect. The position of the null initially moves downwards, before changing direction and reaching the surface after approximately a third of the simulation time. After this point, magnetic field diffuses outwards rapidly. In terms of surface flux, this regime is closer to the 1D case than any other two-step profile we test with the connected active region, though still not a good match.

The middle panel shows field lines from the simulation with the MLT profile and STABLE initial condition. Because of the strong diffusivity in the bulk of the domain, the field spreads out in the convection zone and diffuses radially outwards due to the reduced diffusivity at the surface. This leads to a surface evolution that matches the 1D case very closely.

We now assess the effect of the upper boundary condition on the surface evolution. For this test, we prescribe a constant diffusivity of η = 1 independent of depth. A snapshot of the simulation is shown in the bottom row of Fig. 6. The left-hand panel shows magnetic field lines where the upper boundary condition is potential, and the middle panel shows the field lines where the boundary condition is radial. Qualitative differences are small, but we see in the right-hand panel that there is a little too much magnetic flux at the surface in the radial case, compared to the 1D surface model. Conversely, the potential case matches the 1D evolution closely. If β = 1, we introduce −∂B_x/∂z into Eq. (3) at the surface which is not present in the radial case. Hence the difference between the two regimes is only situated in the upper quarter of the domain. The enforcement of a radial field at the surface boundary also means that the field lines interact with the periodic boundary later because they are strictly vertical, as opposed to the potential case where cancellation can occur more readily. Since the diffusivity is high throughout the domain, field lines can move freely, and the majority of the flux is diffused out of the convection zone by the time we reach t = 0.1.

Although the choice of radial or potential-field boundary condition can slightly change the amount of magnetic flux at the surface, the differences are only small, and starker differences arise when we prescribe a more realistic multi-step diffusion profile in place of the constant diffusivity, or change the connectivity of the active region. Further tests show that the small improvement attained by changing boundary condition is the same regardless of the choice of diffusion profile. Interestingly, the 2D model in the constant case provides a good match to the 1D model and explains in part why the MLT profile performs best out of the multi-step profiles we tested: the strong diffusivity allows the magnetic field to diffuse outwards in both cases, the only difference being that the field lines remain attached to the toroidal field in the MLT case due to a weak base diffusion.

The periodic boundary conditions in x can be interpreted as the presence of neighbouring active regions. To check the influence of this inter-region spacing, we tried increasing the width of the domain. This results in more flux present at the surface because it takes longer to diffuse to the boundary and cancel. However, the results above hold qualitatively, and in any case we cannot choose the locations of active region emergence when simulating the evolution of observed BMRs, so varying the width of the domain does not give us significantly deeper insight.

3. Effect of diffusivity for a 3D decaying active region

We return to the 3D dynamo model KD3 to test whether the results found in Sect. 2 hold qualitatively here as well. We emerge a single region at 10° latitude with flux 1 × 10²² Mx and a tilt angle of 30°. Once the region has emerged after 25 days, the velocity perturbation is turned off. A snapshot of the system is taken on that day, and all subsequent experiments are run from time of emergence. This best reflects the scenario modelled in the simplified 2D diffusion model in Sect. 2, and is the same setup as described in Figs. 1–3 but starting from a different time.

Differential rotation takes the form of Charbonneau et al. (1999):

$$\begin{array}{cc}\hfill \mathrm{\Omega }(r,\theta )=& {\mathrm{\Omega }}_{\mathrm{c}}+\frac{1}{2}[1+\mathrm{erf}\left(\frac{r-R_0}{{\mathrm{\Delta }}_0}\right)][{\mathrm{\Omega }}_{\mathrm{E}}-{\mathrm{\Omega }}_{\mathrm{c}}\hfill \\ \hfill & +({\mathrm{\Omega }}_{\mathrm{p}}-{\mathrm{\Omega }}_{\mathrm{E}})(C{cos}^2\theta +(1-C){cos}^4\theta )],\hfill \end{array}$$(8)

where Ω_c = 2.71434 × 10⁻⁶ s⁻¹, Ω_E = 2.9531 × 10⁻⁶ s⁻¹, Ω_p = 2.07345 × 10⁻⁶ s⁻¹, C = 0.483, R₀ = 0.7 R_⊙ and Δ₀ = 0.025 R_⊙.

For meridional flow we first define the following stream function (Yeates & Muñoz-Jaramillo 2013):

$$\begin{array}{cc}\hfill \mathrm{\Psi }(r,\theta )=& \frac{-v_0(r-R_{\mathrm{p}})}{7.633rsin\theta }sin\left(\pi \frac{r-R_{\mathrm{p}}}{R_{\odot }-R_{\mathrm{p}}}\right)\hfill \\ \hfill & \times exp[-{\left(\frac{r-R_1}{\mathrm{\Gamma }}\right)}^2][1-exp(-1.5{\theta }^2)]\hfill \\ \hfill & \times \{1-exp\left[1.8(\theta -\frac{\pi }{2})\right]\},\hfill \end{array}$$(9)

where R_p = 0.62 R_⊙, R₁ = 0.1125 R_⊙, Γ = 3.47 × 10⁸ m and v₀ = 20 m s⁻¹. Then the meridional circulation is given by

$$\begin{array}{c}\hfill {\mathit{v}}_m=\frac{1}{\rho \left(r\right)}\mathrm{\nabla }\times [\mathrm{\Psi }(r,\theta ){\mathit{e}}_{\varphi }],\end{array}$$(10)

where $\rho \left(r\right)={({\displaystyle \frac{R_{\odot }}{r}}-0.95)}^{\frac{3}{2}}$ is the radial density profile.

We use a grid resolution of ${\mathrm{\Delta }}_{\mathrm{r}}=\frac{0.45R_{\odot }}{48}$ in radius, and a variable grid in latitude and longitude (see Yeates & Muñoz-Jaramillo 2013 for details). Here we set the equatorial grid spacing ${\mathrm{\Delta }}_{\varphi }=\frac{2\pi }{384}$. Initial and boundary conditions are the same as used by Yeates & Muñoz-Jaramillo (2013): the bottom boundary condition at r = 0.55 R_⊙ is a perfectly conducting core, meaning ∂(rB_θ)/∂r = ∂(rB_ϕ)/∂r = 0. The upper boundary condition is radial, although we expect from Sect. 2 that changing to a potential-field boundary condition would have a negligible effect on the flux evolution. The initial condition is created by emerging a single BMR from a layer of toroidal field at the base of the convection zone of the form:

$$\begin{array}{c}\hfill \mathit{B}=-\frac{B_0}{2}\left(\frac{cos\theta }{|cos\theta |}\right)[\mathrm{erf}\left(\frac{r-R_7}{{\mathrm{\Delta }}_8}\right)-\mathrm{erf}\left(\frac{r-R_8}{{\mathrm{\Delta }}_8}\right)]{\mathit{e}}_{\varphi },\end{array}$$(11)

with R₇ = 0.66 R_⊙, R₈ = 0.74 R_⊙, Δ₈ = 0.018 R_⊙ and B₀ = 2.5 × 10³ G. At the surface, we do not prescribe any initial magnetic field, aside from that of the emerged BMR. Diffusivity is now given by η_s η(r/R_⊙) using Eq. (6), where η_s = 6 × 10¹² cm² s⁻¹. We run the model using a different diffusion profile from Sect. 2 each time. The resulting unsigned surface flux and polar flux profiles are shown in Fig. 7.

Fig. 7. Top: unsigned surface flux from 3D simulations of a single active region, using the KD3 diffusion profile (orange), quenching profile (purple) and MLT profile (yellow). The equivalent 2D SFT flux is shown in the dotted blue (without decay) and black (with decay) curves. Bottom: northern polar flux (solid and dotted curves) and southern polar flux (dashed and dash-dotted curves) from the same simulations.

As we found in the 2D model (Sect. 2), the KD3 profile (orange) restricts cancellation, due to the weak diffusivity keeping field lines connected to the toroidal field. This results in a vast excess of flux at the surface. Qualitatively, the other profiles also exhibit the same behaviour as in the 2D model. The MLT profile (yellow) provides a more rapid decay of flux due to the increased diffusivity in the convection zone allowing for disconnection, and the quenching profile (purple) lies somewhere between the other two. By the end of the simulation, there is a similar amount of surface flux in the KD3 simulations as in the SFT simulations, as shown by the dotted blue and black curves. These curves represent the models without the exponential decay term (see Whitbread et al. 2017), and with a decay parameter of τ = 10 years, respectively. Whilst the decay term in the SFT model makes only a very small difference in the total unsigned surface flux, its impact at the poles is more evident, acting as a sink for the polar flux which is not otherwise possible in the SFT model (Baumann et al. 2006). Although the peak strength of the northern polar field is weaker in the MLT case than the SFT model, it occurs at a similar time and the shape of the profile is close to that of the SFT model when exponential decay is included. In summary, these experiments with the decay of a single active region suggest that increasing the diffusivity in the bulk of the convection zone can improve the realism of long-term surface flux evolution compared to the original KD3 model.

4. Effect of diffusivity on a 3D full-cycle simulation

Yeates & Muñoz-Jaramillo (2013) demonstrated a simulation of a full solar cycle using BMR data from Solar Cycle 23. However, this was not systematically calibrated to observations. It can be seen in the top panel of Fig. 8 (or equivalently Fig. 12 of Yeates & Muñoz-Jaramillo 2013) that the magnetic field is too strong and poleward surges are too slow compared to the optimal butterfly diagram found by Whitbread et al. (2017), shown in the lower panel of Fig. 8, which was calibrated against observations. The active regions across the full solar cycle behave similarly to the individual region in Fig. 2. We repeat this 3D simulation of Cycle 23 but replace the original diffusion profile (orange curve in Fig. 5) with the quenching and mixing-length theory profiles (purple and yellow curves in Fig. 5 respectively).

Fig. 8. Top: simulation of Cycle 23 from Yeates & Muñoz-Jaramillo (2013). Bottom: optimal butterfly diagram of Cycle 23 from Whitbread et al. (2017).

Equation (11) again defines the initial toroidal field, but now we try B₀ = 250 G. An initial dipolar field is given by

$$\begin{array}{c}\hfill \mathit{B}=\mathrm{\nabla }\times \left(A_{\varphi }{\mathit{e}}_{\varphi }\right),\end{array}$$(12)

where

$$\begin{array}{c}\hfill A_{\varphi }=B_d\frac{sin\theta }{r^3}\left(\frac{r-0.7R_{\odot }}{0.3R_{\odot }}\right),\end{array}$$(13)

and A_ϕ = 0 for r <  0.7 R_⊙ (Jouve et al. 2008). The field strength is set as B_d = −0.008 B₀.

We run the simulation for 5000 days, using observed BMRs of Cycle 23 from NSO Kitt Peak as input data (Yeates et al. 2007), as in Yeates & Muñoz-Jaramillo (2013). The unsigned surface flux and signed polar flux for the simulation of Yeates & Muñoz-Jaramillo (2013) are shown by the orange curves (top and bottom respectively) in Fig. 9. The purple profiles in this plot correspond to the simulation where the quenching diffusivity profile has been used. If all parameters other than the diffusivity profile are fixed, it is evident that not enough magnetic flux reaches the surface, and the polar field is barely able to reverse.

Fig. 9. Top: unsigned surface flux from 3D simulations of Cycle 23 using the KD3 diffusion profile (orange), quenching profile (purple) and MLT profile (yellow). Bottom: northern polar flux (solid line) and southern polar flux (dashed line) from the same simulations.

To combat this, we increase the strength of the initial toroidal field by an order of magnitude. This provides a stronger source from which active regions can develop, thereby increasing the amount of flux at the photosphere. This is demonstrated by the purple curve in the top panel of Fig. 10. Here, the total surface flux peaks earlier than the original simulation. In the bottom panel, we see that the polar field reverses at a similar time to the original case, albeit with a reduced strength throughout the simulation. Nevertheless, the toroidal field appears to be strong enough to produce more regions as a subsequent cycle (top panel of Fig. 11) if we were to continue the simulation. The bottom panel of Fig. 11 shows the surface butterfly diagram of the same simulation. While it is suboptimal, it displays observable features of the solar cycle and a more realistic distribution and transport of magnetic flux than before. A future task is to calibrate other parameters in the model against observations while keeping the quenching profile fixed.

Fig. 10. Top: unsigned surface flux from 3D simulations of Cycle 23 using the KD3 diffusion profile (orange), quenching profile (purple) and MLT profile (yellow), but where the initial toroidal field has been strengthened by one and two orders of magnitude for the latter two respectively. Bottom: northern polar flux (solid line) and southern polar flux (dashed line) from the same simulations.

Fig. 11. Top: toroidal field at the base of the convection zone from a 3D simulation of Cycle 23 using the quenching profile and a strengthened initial toroidal field. Bottom: radial magnetic field at the surface from the same simulation.

Ideally we would like to be able to simulate Cycle 23 using the diffusion profile derived from mixing-length theory, because the enhanced diffusivity acts as a means of disconnecting the emerged regions from the toroidal field, and this profile gave the closest match to the surface-only evolution in Sects. 2 and 3. Figure 9 shows that even less flux emerges at the surface in this case, because the diffusion in the convection zone is too strong and kills off the majority of rising flux tubes and returning poloidal flux. Even when the initial toroidal field is increased by an order of magnitude, it rapidly diffuses and so no regions are able to emerge after a few years.

When the toroidal field is increased by another order of magnitude, the flux still decays too rapidly, as shown by the yellow curve in Fig. 10. However, we now observe polar field reversal, although very early in the cycle, and the bottom panel of Fig. 12 shows that the surface evolution during the first few years of the cycle appears to be sun-like. The top panel of Fig. 12 shows that no new toroidal field is created. This occurs in all simulations when the MLT diffusivity profile is used and is one reason why dynamos have thus far been unable to accommodate the diffusion profile derived from mixing-length theory.

Fig. 12. Top: toroidal field at the base of the convection zone from a 3D simulation of Cycle 23 using the MLT profile and a strengthened initial toroidal field. Bottom: radial magnetic field at the surface from the same simulation.

Scaling the MLT profile by a factor of 0.5 allows significantly more flux to emerge at the surface, but it is still not enough on its own to sustain the dynamo. However, if we also shift the location of the low-diffusivity step in the MLT profile up so that the toroidal field is stored in a region of low diffusion (i.e. set R₁ = 0.74 and Δ₁ = 0.024), we find that the field survives for longer and more flux can reach the surface. However, although more new toroidal field starts to appear at the base of the convection zone for the next cycle, it is still too weak, and the polar field at the surface still reverses too early. In summary, increasing the diffusivity to the level required for a realistic surface evolution is not on its own sustainable in a full-cycle simulation, because the high diffusivity removes too much flux from the system.

5. Conclusions

Our main conclusion is that 3D kinematic dynamo models cannot produce a realistic evolution of magnetic flux on the solar surface if active regions are allowed to remain connected to the base of the convection zone. Within the framework of the KD3 model, where active regions are formed self-consistently through imposed velocity perturbations, the only way to achieve disconnection is through enhanced turbulent diffusivity in the convection zone. Whilst such an enhanced diffusivity is compatible with estimates from mixing-length theory, flux transport dynamo models have been unable to function with such a high diffusivity (Muñoz-Jaramillo et al. 2011). In this paper we have demonstrated that, indeed, a full-cycle simulation with KD3 is not possible with such strong diffusivity, despite the fact that it leads to a realistic surface flux evolution when simulating the decay of a single active region.

A possible resolution to this problem is suggested by simulations where the active region is initially disconnected from the base of the convection zone, as in the 3D dynamo model STABLE of Miesch & Dikpati (2014). With magnetic field lines no longer anchored to the base of the convection zone, diffusion is much more effective. We have shown clearly in our 2D simplified model (Sect. 2) that this leads to a better fit with the surface evolution, for a given diffusivity profile. In the bottom panel of Fig. 13, we demonstrate the effect of disconnecting regions in KD3, using the quenching diffusion profile as an illustration (purple curves). For the interior magnetic field we calculate a potential field extrapolation inward from the given surface B_r with a perfectly conducting boundary condition at a fixed depth. We find that the flux decays faster than in the SFT case, but that the evolution is dependent on the depth of the potential field. If the region is shallow (solid curve), it decays very quickly, but if we increase the depth to 0.8 R_⊙ (dash-dotted curve) and then 0.65 R_⊙ (dashed curve), we observe an increasingly better fit to the 2D surface evolution, although it is still by no means perfect. With further study, it may be possible to produce an excellent fit to the surface model using disconnected regions. The top row of Fig. 13 shows two examples of disconnected regions from KD3 with depths 0.65 R_⊙ (left) and 0.95 R_⊙ (right). The shallower depth forces the outermost field lines to be pressed close to the surface (cf. Fig. 1).

Fig. 13. Top: three-dimensional images of a disconnected active region in KD3 with depth 0.65 R⊙ (left) and 0.95 R⊙ (right). The radial magnetic field is shown at the transparent surface and field lines below the surface. Bottom: unsigned surface flux from 3D simulations of a single disconnected active region with depth 0.65 R⊙ (dashed), 0.8 R⊙ (dash-dotted) and 0.95 R⊙ (solid), using the quenching profile. The equivalent 2D SFT flux is shown in the dotted blue (without decay) and black (with decay) curves.

Several observational facts at the surface point to the likelihood that active regions become rapidly disconnected from their roots after emergence, on a timescale of days (Fan 2009). For example, if they remained connected we might expect to see a relaxation of tilt angle towards the east-west line once emergence (and the resulting Coriolis force) ceased to operate, as well as a continued separation of the two polarities reflecting conditions at the roots of the emerged flux tube. Indeed, Fan et al. (1994) remarked that the success of SFT models in reproducing the surface evolution is itself evidence that active regions are no longer dynamically constrained from below. The actual process of disconnection is less understood, although Schrijver & Title (1999) argued that subsurface reconnection on the required timescale will be a natural consequence of convective motions. Fan et al. (1994) proposed a mechanism of “dynamical disconnection”, based on the idea that the magnetic field in the subsurface legs of the active region could be weakened to sub-equipartition values as it tries to establish hydrostatic equilibrium, thus enabling reconnection. This was confirmed in a 1D calculation by Schüssler & Rempel (2005).

Let us end with some broader remarks. In this paper, we have treated the SFT model as the “ground truth”, but in reality we must remember that 2D SFT models have their own limitations. The resulting errors are likely too small to change our broad conclusions in this paper, but may be non-negligible. For example, Fig. 7 shows how adding an exponential decay term to the SFT model improves the qualitative match with the 3D model when viewed over solar-cycle timescales. This term is a crude parametrization of radial diffusion, which is missing in the basic SFT model (Baumann et al. 2006) but included consistently in the 3D model. Our results suggest that, for the diffusivities used in typical flux-transport dynamo simulations, this term will have a non-negligible effect on the surface evolution.

At this point, the important question remains as to whether a flux-conserving 3D B–L model like KD3 could match the observed surface flux evolution on the Sun. From our results, we expect that the existing dynamo codes with non-local active region deposition, like STABLE, are better able to match the surface evolution than the existing KD3 code. However, KD3 nevertheless has the advantage of satisfying the magnetic induction equation during the emergence of active regions, allowing the critical flux budget to be correctly accounted for. In future, it would therefore be desirable to develop an active region emergence model that satisfies the induction equation while at the same time achieving rapid disconnection of active regions, perhaps through controlled reconnection. Exploring the full parameter space of such a model, including the meridional flow and turbulent pumping profiles, would be a significant computational task. Nevertheless, the work in this paper, along with restrictions on the surface parameters from Whitbread et al. (2017), does add some valuable constraints to this optimization problem.

Acknowledgments

TW thanks Durham University Department of Mathematical Sciences for funding his PhD studentship. ARY thanks STFC for financial support through grant ST/N000781/1. This research was supported by the NASA Living With a Star Grant NNH15ZDA001N.

References

All Figures

	Fig. 1. Three-dimensional image of an emerged active region in KD3. Magnetic field lines are connected to the toroidal field at the base of the convection zone and the radial magnetic field is shown at the transparent surface.
In the text

	Fig. 2. Top: longitude-averaged evolution of Br for a single BMR in a 2D SFT model. Bottom: surface component of the 3D dynamo model showing the equivalent evolution of the same BMR.
In the text

	Fig. 3. Top: comparison of unsigned surface flux from the 2D SFT simulation (blue) and 3D dynamo simulation (orange). Bottom: comparison of northern (solid and dotted lines) and southern (dashed and dash-dotted) polar flux from the same two simulations, where polar flux is defined as the flux polewards of 70° latitude.
In the text

	Fig. 4. Two initial conditions used in this paper: an active region connected to the toroidal field (left) and an active region disconnected from the toroidal field (right).
In the text

	Fig. 5. Normalised multi-step diffusion profiles used in this paper, against a log-scale. The solid orange curve is from KD3, the dashed purple curve is the profile which takes into account diffusivity quenching, and the dotted yellow curve is derived from mixing-length theory.
In the text

Fig. 6. Snapshots of magnetic field lines from simulations using the KD3 (top row), quenching (second row), MLT (third row) and constant (bottom row) diffusion profiles. For the first three profiles, the left-hand column shows the case where the region is connected to the base, and the middle column shows the simulation with a disconnected initial condition. The black dashed line is the top of the domain, above which is shown a potential-field extrapolation. Right-hand column: evolution of magnetic flux at the surface, compared to a 1D surface model (green dashed line). The bottom row (constant η) compares the effects of potential (blue) and radial (red) boundary conditions but using the same connected initial condition.

In the text

Fig. 7. Top: unsigned surface flux from 3D simulations of a single active region, using the KD3 diffusion profile (orange), quenching profile (purple) and MLT profile (yellow). The equivalent 2D SFT flux is shown in the dotted blue (without decay) and black (with decay) curves. Bottom: northern polar flux (solid and dotted curves) and southern polar flux (dashed and dash-dotted curves) from the same simulations.

In the text

	Fig. 8. Top: simulation of Cycle 23 from Yeates & Muñoz-Jaramillo (2013). Bottom: optimal butterfly diagram of Cycle 23 from Whitbread et al. (2017).
In the text

	Fig. 9. Top: unsigned surface flux from 3D simulations of Cycle 23 using the KD3 diffusion profile (orange), quenching profile (purple) and MLT profile (yellow). Bottom: northern polar flux (solid line) and southern polar flux (dashed line) from the same simulations.
In the text

	Fig. 10. Top: unsigned surface flux from 3D simulations of Cycle 23 using the KD3 diffusion profile (orange), quenching profile (purple) and MLT profile (yellow), but where the initial toroidal field has been strengthened by one and two orders of magnitude for the latter two respectively. Bottom: northern polar flux (solid line) and southern polar flux (dashed line) from the same simulations.
In the text

	Fig. 11. Top: toroidal field at the base of the convection zone from a 3D simulation of Cycle 23 using the quenching profile and a strengthened initial toroidal field. Bottom: radial magnetic field at the surface from the same simulation.
In the text

	Fig. 12. Top: toroidal field at the base of the convection zone from a 3D simulation of Cycle 23 using the MLT profile and a strengthened initial toroidal field. Bottom: radial magnetic field at the surface from the same simulation.
In the text

Fig. 13. Top: three-dimensional images of a disconnected active region in KD3 with depth 0.65 R⊙ (left) and 0.95 R⊙ (right). The radial magnetic field is shown at the transparent surface and field lines below the surface. Bottom: unsigned surface flux from 3D simulations of a single disconnected active region with depth 0.65 R⊙ (dashed), 0.8 R⊙ (dash-dotted) and 0.95 R⊙ (solid), using the quenching profile. The equivalent 2D SFT flux is shown in the dotted blue (without decay) and black (with decay) curves.

In the text