Velocity Amplitudes in Global Convection Simulations: The Role of the Prandtl Number and NearSurface Driving
Abstract
Several lines of evidence suggest that the velocity amplitude in global simulations of solar convection, , may be systematically overestimated. Motivated by these recent results, we explore the factors that determine and we consider how these might scale to solar parameter regimes. To this end, we decrease the thermal diffusivity along two paths in parameter space. If the kinematic viscosity is decreased proportionally with (fixing the Prandtl number ), we find that increases but asymptotes toward a constant value, as found by Featherstone and Hindman (2016). However, if is held fixed while decreasing (increasing ), we find that systematically decreases. We attribute this to an enhancement of the thermal content of downflow plumes, which allows them to carry the solar luminosity with slower flow speeds. We contrast this with the case of RayleighBénard convection which is not subject to this luminosity constraint. This dramatic difference in behavior for the two paths in parameter space (fixed or fixed ) persists whether the heat transport by unresolved, nearsurface convection is modeled as a thermal conduction or as a fixed flux. The results suggest that if solar convection can operate in a high regime, then this might effectively limit the velocity amplitude. Smallscale magnetism is a possible source of enhanced viscosity that may serve to achieve this high regime.
keywords:
Sun: interior, convection[cor]Corresponding author
1 Introduction
Observations clearly indicate that the surface of the Sun is blanketed by roughly a million convection cells, each with a characteristic size of about 1.4 Mm and a typical lifetime of about 1015 min. This is the wellknown phenomena of solar granulation and it represents the dominant scale of convection in the solar photosphere. Granulation is well understood theoretically and modern numerical simulations are able to reproduce observations with striking fidelity (Nordlund et al., 2009). In fact, modeling solar granulation is arguably one of the most notable success stories in all of computational astrophysics.
However, there is far more to solar convection than solar granulation. Granulation alone cannot account for the differential rotation and meridional circulation apparent from surface observations and helioseismic inversions. Furthermore, it is these mean flows, together with the deeper, largerscale convective motions that must maintain them, that are thought to be responsible for establishing the solar activity cycle (Miesch, 2005; Charbonneau, 2010).
There are three recent research developments that suggest our current understanding of deep solar convection may be incomplete. The first comes from the same numerical simulations of solar surface convection that do so well at capturing granulation. When these are extended to larger horizontal and vertical scales, they tend to overestimate the power at low spectral wavenumbers when compared to photospheric observations (Lord et al., 2014). For the largest horizontal scales in the simulations, 196 Mm, the characteristic amplitude of convective velocities is at least three times larger than that inferred from correlation tracking of photospheric features, with the trend suggesting even larger discrepancies at larger scales. Furthermore, the peak vertical velocities at the bottom of the computational domain, 49 Mm below the solar surface, can exceed several hundred m s. Though these surface convection simulations neglect rotation, such large convective velocities suggest a weak rotational influence.
To demonstrate this, we can quantify the rotational influence by means of the Rossby number where is some measure of the rotation rate and and are typical convective velocity and length scales. If we take rad s 200 m s, and 28 Mm, where is the density scale height, we obtain . This appears to be at odds with the differential rotation of the Sun inferred from photospheric observations and helioseismic inversions, which reveal a 30% decrease in angular velocity from equator to pole. Such a substantial differential rotation implies efficient angular momentum transport by the convective Reynolds stress (Miesch, 2005). This in turn requires the strong Coriolisinduced velocity correlations that are generally associated with .
This expectation is borne out by global simulations of solar and stellar convection, which show a dramatic change in the differential rotation profile near (Guerrero et al., 2013; Gastine et al., 2014; Fan and Fang, 2014; Käpylä et al., 2014; Featherstone and Miesch, 2015; Karak et al., 2015). For , differential rotation profiles are generally solarlike in the sense that the gradient is equatorward, decreasing from equator to pole (though the contours are often more cylindricallyaligned than in the Sun). For , differential rotation profiles are antisolar, with a poleward gradient. This is attributed to the tendency for convective motions to conserve their angular momentum locally when the rotational influence is weak.
Global simulations of solar convection often fall in the antisolar rotation regime, suggesting that they are overestimating the Rossby number of the deep solar convection zone. This contradiction is consistent with the velocity amplitudes of surface convection simulations described above but represents a second, independent recent research development that calls into question our current understanding of deep solar convection. Current global simulations of solar convection often artificially reduce the Rossby number by increasing the rotation rate (Augustson et al., 2015; Karak et al., 2015) or by decreasing the convective luminosity, which in turn suppresses (Racine et al., 2011; Guerrero et al., 2013; Fan and Fang, 2014; Hotta et al., 2015a).
Recent results from local helioseismology also suggest that numerical and theoretical models may be overestimating the velocity amplitudes of deep solar convection, though the evidence is mixed. This is the third development referred to above. Hanasoge et al. (2010, 2012) searched for signatures of giant cells in helioseismic timedistance measurements and did not find a significant signal. They placed an upper limit of less than 3 m s for each spherical harmonic degree at for convective motions with coherence times longer than 24 hours and , where is the solar radius. When the coherence time was increased to 96 hours, they inferred amplitudes of less than 1 m s for each . When summed over all spherical harmonic modes , this corresponds to an upper limit on of less than 10 m s. This upper limit is more than an order of magnitude lower than the convective velocities typically found in global convection simulations and predicted from mixinglength theory. Furthermore, it is hard to reconcile with helioseismic measurements of the differential rotation and meridional circulation (Miesch et al., 2012). More recent helioseismic estimates of the convective flow speed using ringdiagram analysis by Greer et al. (2015) suggest values that are closer to the theoretical models. They find 120 m s at 0.960.97 . For a recent review of these and other results, see Hanasoge et al. (2016).
Though each of these three recent developments is subject to its own uncertainties, they collectively suggest that there may be something fundamental about deep solar convection that we are missing. In short, we may state the convective conundrum as follows: The convective velocities required to transport the solar luminosity in global models of solar convection appear to be systematically larger than those required to maintain the solar differential rotation and those inferred from solar observations.
Even more to the point of this paper, the problem seems to get worse as the dissipation is decreased. Due to computational constraints, the values of the kinematic viscosity and the thermal diffusivity that are currently used in global convection simulations are orders of magnitude larger than the corresponding values of the actual solar plasma. One would hope that the conundrum might resolve itself if and could be drastically reduced, making the parameter regime closer to the Sun. However, simulations have instead found that lowering and can exacerbate the problem by further increasing the effective Rossby number (Guerrero et al., 2013; Gastine et al., 2014; Featherstone and Miesch, 2015; Karak et al., 2015).
Motivated by this conundrum, we have performed idealized numerical experiments to investigate the nature of the convective driving in global simulations. We focus in particular on significance of the Prandtl number, which is the ratio of viscous to thermal diffusion; . We demonstrate that a high value of helps to limit as is decreased and we argue that this may help to resolve the convective conundrum. We also demonstrate that this conclusion is insensitive to the details of how the heat transport by unresolved nearsurface convection is parameterized.
In section 2 we describe the Rayleigh code, which we use to solve the anelastic fluid equations in a spherical shell, as well as the setup of the numerical experiments. For these fundamental studies we neglect rotation. We present the results of our numerical experiments in sections 3 and 4, which focus on the role of the Prandtl number and the upper boundary layer respectively. We close with a discussion and summary of our results and their implications in sections 5 and 6.
2 Numerical Experiments
2.1 Numerical Model
This study is based around a series of 3D, nonlinear convection models run with the Rayleigh convection code. Rayleigh is a modern, parallel code that has been designed to simulate thermal convection in spherical geometries at high resolution; see Featherstone and Hindman (2016, hereafter FH16) and a forthcoming code paper for further details.
Rayleigh employs a pseudospectral approach, meaning that derivatives along each axis are calculated in the relevant spectral configuration, whereas nonlinear terms are calculated after first transforming to physical space. We represent the horizontal variation of all variables along spherical surfaces using spherical harmonics where is the spherical harmonic degree, and is the azimuthal order. In the radial direction, we expand all variables using Chebyshev polynomials , where is the degree of the Chebyshev polynomial. Timeintegration is carried out in the spectral configuration using a hybrid implicitexplicit approach. A CrankNicolson method is used for linear terms, and an AdamsBashforth approach is employed for the nonlinear terms. Both components of the timestepping are 2ndorder in time.
The numerical algorithm for Rayleigh is similar to the wellknown Anelastic Spherical Harmonic (ASH) code (Clune et al., 1999; Brun et al., 2004). However, Rayleigh has been designed to be a bit more flexible (with a Cartesian option) and more scalable (somewhat higher parallel efficiency) than ASH. It is also somewhat more accessible for both students and the community; the plan is to distribute it publicly by means of the Computational Infrastructure for Geophysics (CIG; see http://geodynamics.org/cig).
Our study is concerned with convection as it manifests deep within stellar interiors, far removed from the photospheric surface where radiative processes may contribute considerably to the energetics of the convection. In such a region of the star, plasma motions are subsonic and perturbations to thermodynamic variables are small compared to their mean, horizontally averaged values. Under such conditions, the anelastic approximation, which we employ, provides a convenient means of describing the thermodynamic background state of the system (Gough, 1969; Gilman and Glatzmaier, 1981).
Under this approximation, thermodynamic variables are linearized about a spherically symmetric reference state with density , pressure , temperature , and specific entropy . Here we use an adiabatic reference state so . Fluctuations about this state are denoted as , , , and . A further consequence of the anelastic approximation is that the mass flux is solenoidal, reducing the continuity equation to
(1) 
where is the velocity vector expressed in spherical coordinates: . The lack of any time derivative in Equation (1) means that sound waves are naturally filtered out as a consequence of this approach. The divergencefree constraint for the mass flux is enforced by projecting onto poloidal and toroidal streamfunctions ( and respectively), such that
(2) 
The unit vector in the radial direction is indicated by . The momentum equation is given by
(3) 
where is the gravitational acceleration. The viscous stress tensor is given by
(4) 
where is the strain rate tensor, the kinematic viscosity is denoted by , and is the Kronecker delta. Our thermal energy equation is given by
(5) 
where the thermal diffusivity denoted by . Both and are assumed to be constant, independent of , , , and . Sources of internal heating and cooling are encapsulated by (see Sec. 2.2). A linearized equation of state closes our set of equations and is given by
(6) 
assuming the ideal gas law
(7) 
The specific heat at constant pressure is denoted by , is the gas constant, and is the adiabatic index of the gas.
2.2 The Numerical Experiments
We have constructed a series of fourteen stellar convection zone models designed to explore how the convective kinetic energy scales as we decrease the dissipation, approaching solar parameter regimes as discussed in Sec. 1. We focus in particular on how this scaling depends on the Prandtl number , and on the prescription for cooling at the upper boundary. These simulations are summarized in Table 1.
Each of our models is constructed using a polytropic background state following Jones et al. (2011). This approach has the advantage that the thermodynamic background may be specified analytically, making it easily reproducible. All the simulations reported here have the same background state, which is completely specified by seven numbers: the inner radius of the shell cm (0.718 ), the outer radius of the shell cm (0.946), the polytropic index , the number of density scale heights occurring within the shell , the mass interior to the shell g, the density at the inner boundary g cm, and erg g K. We note that when constructed using these parameters, the thermodynamic background state closely resembles that of the Sun (FH16).
For all simulations, we have adopted impenetrable and stress free boundary conditions such that
(8) 
For the C series of simulations listed in Table 1, the radial entropy gradient is forced to vanish at the lower boundary of the convection zone, and the entropy perturbations are forced to vanish at the upper boundary:
(9) 
Thus, there is no diffusive entropy flux across the lower boundary. Heat is deposited into this system by instead, which drops to zero at the upper boundary. In all simulations, we adopt a functional form for that depends only on the background pressure gradient such that
(10) 
The normalization constant is chosen so that
(11) 
where is the solar luminosity. We then define a corresponding radial energy flux as follows:
(12) 
where . We give this flux the subscript because we identify it with the flux due to radiative diffusion. It is defined such that , and . Thus, radiative diffusion carries the entire solar luminosity in to the computational domain through the bottom boundary and no flux through the top boundary. This implies a convergence of that establishes and sustains the convection.
To illustrate how this occurs, we consider the other components of the radial heat flux, which include the convective enthalpy flux , the kinetic energy flux , and the conductive flux , defined as follows:
(13)  
(14)  
(15) 
Angular brackets denote averages over horizontal surfaces ( and ). In a steady state, these fluxes must balance such that .
Each simulation is initialized using a small random thermal fluctuation. As time proceeds, the convergence of the radiative heat flux heats the convection zone (CZ), raising the adiabat and establishing a negative near the surface, as illustrated in the right panel of Fig. 1. This superadibatic entropy gradient excites convective motions that grow to dominate the heat transport through the mid CZ by means of the resolved convective fluxes and . The resulting balance is shown in Fig. 1 for Case C1 as an example. Each simulation was evolved until the kinetic energy reached a statistically steady (equilibrated) state, and then further evolved for at least one thermal diffusion time (some were evolved up to ten thermal diffusion times).
For the C series of simulations, heat exits the upper boundary via the conductive flux, . There, the steepness of the equilibrated entropy gradient is entirely dependent upon the thermodynamic background state, the value of the thermal diffusivity , and the chosen luminosity. Specifically,
(16) 
This leads to a variation among simulations in the entropy difference across the CZ, . Equilibrated values of for all simulations are included in Table 1.
For the G series of simulations, eq. (16) no longer holds. There we have imposed another flux component that contributes to the heat transport through the outer surface (see Table 1). This is intended to mimic heat transport by smallscale surface convection as discussed in §4. We have also modified the upper thermal boundary condition in the G series, fixing the entropy gradient, , as opposed to the entropy , as expressed in eq. (9). Apart from these differences, the series was initiated and evolved in a similar manner as the series. The difference between the different members of the C and G series are the values of and , as indicated in Table 1.
Case  Surface  

( cm s)  ( cm s)  Flux  erg g K  
C1  8  8  1  4620  
C2  6  8  1.33  5560  
C3  4  8  2  7060  
C4  2  8  4  10300  
C5  6  6  1  5450  
C6  4  4  1  6780  
C7  2  2  1  8870  
G1  8  8  1  3040  
G2  6  8  1.33  2690  
G3  4  8  2  2422  
G4  2  8  4  2033  
G5  6  6  1  2704  
G6  4  4  1  2312  
G7  2  2  1  1900 
3 Results: The Role of the Prandtl Number
The Prandtl Number () is defined as (Sec. 2.1). FH16 investigated the scaling of the kinetic energy with in global, nonrotating convection simulations for the special case in which the value of is fixed at unity. They found that the characteristic velocity amplitude first increased slightly as was decreased but then asymptoted toward a constant value in the limit . Here we consider the case in which is held fixed as is decreased (implying an increasing value of ) and we find that in this case systematically decreases.
In figure 1 (left) we show the flux balance across the convection zone for an illustrative simulation in the C series, case C1. The figure shows that is the main transporter of the solar luminosity at low values of (see Sec. 2.2). Near the center of the convection zone () transports nearly all the solar luminosity to the sun’s surface, while at the boundaries ( and ) transports none of the solar luminosity. takes over as the main carrier of the solar luminosity as .
In the right panel of figure 1 we look at the specific entropy over for various values. As the value of decreases (blue being case C1 and tiel being case C4) we see an increase in the change of the specific entropy from the bottom to the top of the convection zone. As noted in Sec. 2.2, this change occurs because a lower value causes an increase in the entropy gradient between the center and top of the convection zone (see eq. 16).
Figures 2 and 3 compare the structure of the convection in two simulations with the same value of but different values of (Cases C1 and C4). A lower value increases both the Rayleigh number and the Prandtl number, causing the convection to become more intricate and apparently more turbulent (Fig. 2). Both the characteristic size of the convection cells and the characteristic width of downflow lanes is smaller for Case C4 (right panel) than for Case C1 (left panel). The latter can be attributed to the thinner thermal boundary layer (see the right panel of Fig. 1 and Sec. 5).
Similar results are observed in the snapshots of the specific entropy in the mid CZ for cases C1 and C4 (Fig. 3). The smaller value of yields a more intricate, smallerscale structure with larger thermal fluctuations. In the right panel of Fig. 3 smaller and sharper features can be seen than in the left panel. This is particularly noticeable for the cold spots (blue) which represent small, cool plumes of fluid that sink down from the surface layers.
Also notable in Fig. 3 is a pronounced dipole convective mode, visible here as a pattern with longitudinal wavenumber . This is a common feature for convection in nonrotating spherical geometries, related to the linearly preferred convective modes (Chandrasekhar, 1961). It is analogous to the box mode or “wind” that commonly arises in studies of convection in Cartesian geometry. Its presence is further promoted by the fixedheat flux lower boundary condition, which favors convective modes with long horizontal wavelengths (Chapman and Proctor, 1980; Depassier and Spiegel, 1981). The dipolar mode is most visible in temperature or entropy plots but it can be seen in radial velocity plots as well (e.g. Fig. 2). The primary effect of the dipole mode is to sweep smallscale, plumelike structures out of the dipolar upflow region, thereby causing them to cluster in the region of convergence associated with the dipolar downflow region. This effect is not observed in rotationally constrained convection because the presence of rotation tends to favor convective modes with higher azimuthal wavenumber (high banana cells).
However, due to the spherical symmetry of the equations, the orientation of this dipole mode should be random; there is no particular reason for it to be nearly perpendicular to the rotation axis as it is in both cases highlighted in Fig. 3 (see also Fig. 4 below). This may be a coincidence; other simulations (not shown here) do exhibit a wide range of orientations for the dipole mode. However, it’s also possible that the numerical grid introduces a slight symmetry breaking that is enough to influence the orientation of the dipole mode. In any case, we do not expect the orientation of the dipole mode to significantly influence our results. The faint horizontal striping pattern in Fig. 3 appears to be an artifact arising from the python plotting program; it is not present in the raw data (before the Molleweide projection).
If is decreased along with , this also leads to a flow that appears more turbulent. This can be seen by comparing the left panel of Fig. 4 to Fig. 2 and the right panel of Fig. 4 to Fig. 3. This greater level of turbulence is usually attributed to an increase of the Rayleigh number, which scales as (FH16).
The apparent increase in the level of turbulence for cases with lower dissipation also suggests a greater convective amplitude, . However, in Fig. 5 we see that this is not necessarily the case. In the left panel we have plotted the root mean square (rms) value of the convective velocity versus 1/ for the C series of simulations. We use this rms velocity as a measure of . The results for fixed agree with those presented by FH16 and show increasing toward a constant value as is decreased. However, the results for fixed show a decrease in with decreasing .
The level of turbulence is typically quantified by the Reynolds number, , where is a characteristic length scale of the flow. If we take to be the depth of the layer, then a decrease of at fixed also implies a decrease in . This conclusion also holds if we use a more meaningful measure of the turbulent length scale , such as the integral scale of the flow or the Taylor microscale. Either of these measures would imply a decrease in with decreasing , further reducing . This can be verified by considering the power spectra shown in the right panel of Fig. 5. As decreases the power spectrum shifts toward higher values of the wave number (). This indicates smaller eddies, which in turn implies a decrease in .
Thus, we arrive at the somewhat surprising conclusion that the flow at low values of and fixed is actually less turbulent (in the sense of smaller ), despite its more intricate, smallerscale structure, and higher Rayleigh number.
4 Heat Transport by SmallScale Surface Convection
Since the entropy fluctuations responsible for the buoyant driving of convection originate mainly in the upper thermal boundary layer, one might expect convective velocity amplitudes to be sensitive to the details of how this region is modeled (see §5). In the Sun, heat transport in the outer 10 Mm of the convection zone () is dominated by relatively smallscale convective motions ranging in size from about 12 Mm (granulation) to about 3035 Mm (supergranulation). Near the solar photosphere, convective transport transitions to radiative heat transport that carries energy away from the solar surface and into interplanetary space. No existing global solar convection simulation has the spatial resolution, temporal resolution, and physical ingredients (radiative transfer, nonideal equation of state, compressibility) to realistically capture the complex dynamics in this outermost region of the Sun. So, we must rely on modeling.
The approach we have taken for the simulations presented in §3 parallels that taken by many previous studies, dating back to the pioneering work of Gilman (1983) and Glatzmaier (1985); see Miesch (2005) and Jones et al. (2011) for further references and discussion. In particular, we modeled the heat flux by unresolved, smallscale convection in the near surface layers as a turbulent thermal diffusion, , that operates on the entropy gradient as opposed to the temperature gradient (see eqs. 5 and 15). This mimics the tendency for efficient convection to mix entropy, establishing a nearly adiabatic stratification.
Though physically justified, the representation of subgridscale (SGS) heat transport as a thermal conduction is an approximation that needs to be evaluated. One consequence of this model is that it ties the SGS heat flux to the background entropy gradient, , which has important implications for the overall dynamics of the convection zone (see §5).
In order to assess these implications further, we have initiated another series of simulations with a different model for the SGS heat flux in the surface layers. This model consists of an imposed radial heat flux, , that is independent of the background stratification and the convective flow field. Its functional form is given by
(17) 
If the parameter is set to unity, then this formulation stipulates that will carry the entire solar luminosity in the region . For , . The width of the transition at is governed by the parameter . Here we use and .
For this series of simulations, which we refer to as the G series (), we also use a different upper thermal boundary condition. Instead of fixing as in the series, we fix the entropy gradient at a value of cm K s. Since we do not zero out , this implies that we we still have a nonzero thermal conduction, , that contributes to the heat flux through the top boundary. However, since is fixed, as .
This brings us to the significance of the parameter . This is an amplitude factor that we adjust to ensure that the total flux through the outer boundary () is equal to the solar flux; . As , more of the flux through the boundary must be carried by , so . However, for the values of considered here, the contribution from is still significant, implying .
This is illustrated in the left panel of Fig. 6, which shows the flux balance for Case G3. Here, transports about 28% of the solar luminosity through the surface, leaving about 72% to be carried by (). The right panel shows the resulting entropy stratification. Here the spread in is smaller than in the C series and shows an opposite trend, decreasing slightly with decreasing instead of increasing. Furthermore, unlike the C series, the entropy profile in the G series will become independent of in the limit .
Remarkably, even with this very different treatment of the upper boundary layer, the results are similar. This is demonstrated for the rms velocity amplitude in the left panel of Fig. 7. Though the velocity amplitudes are systematically higher in the G series of simulations, the trends are similar to the C series. In particular, decreases with decreasing if is held fixed but increases and asympotes toward a constant value if is held fixed. The velocity spectra are also similar, though somewhat steeper in the G series (Fig. 7, right panel). This can likely be attributed to the greater width of the thermal boundary layer for nonzero , which is expected to produce wider plumes (§5).
Thus, we conclude that the results, and in particular the velocity amplitudes, are insensitive to the details of how the upper thermal boundary layer is modeled. This is consistent with the highresolution, global, nonrotating simulations of Hotta et al. (2014). They compared several simulations with different surface cooling and different locations for the upper boundary, and . The latter employed an unprecedented spatial resolution in order to explicitly capture smallscale convective motions near the surface (though even this simulation did not have enough resolution to extend the convective spectrum all the way down to granulation). They found that various measures of the velocity field in the mid convection zone, including the velocity amplitude and spectrum, were similar in the different cases considered and thus insensitive to the detailed structure and dynamics of the upper boundary layer.
5 Discussion
When interpreting our results, we must keep in mind the target parameter regime. As discussed in §4, the conductive flux is a model that is intended to represent SGS convection in the solar surface layers. As such, there is no value of that corresponds to the microphysical (plasma) thermal diffusion in the actual Sun. However, there is a thermal diffusion due to radiative heat transport (modeled here implicitly by ), which is proportional to . Thus, we may loosely take the value of cm s as our target value for the effective in the Sun. Since the values of used here are on the order of  cm s, we would need to decrease by about 56 orders of magnitude to approach the actual thermal diffusion of the Sun. The molecular viscosity, , in the Sun is even smaller, on the order of 1 cm s (Miesch, 2005). So, in order to achieve the true microphysical plasma viscosity of the Sun, we would have to decrease by about 1213 orders of magnitude. So, based on the microphysics of the solar plasma, . However, as we will argue in §6, the effective viscosity of the solar plasma may be much higher, potentially placing the Sun in a high regime.
This justifies our two paths in parameter space. We seek to understand the behavior of the system as we decrease to approach more solarlike parameter regimes. However, we find that this behavior depends on whether or not we simultaneously decrease as well.
Within this context, the aim of this section is to understand and interpret the principle result of Section 3, namely that:

R1: decreases with decreasing if is held fixed.
where is the characteristic velocity amplitude of the convection. If this result holds up then it bodes well for global solar convection simulations, that employ an artificially large value of and that appear to be overestimating (§1).
While we are addressing this principle result, we will also address the other results discussed in sections 3 and 4, which may be summarized as follows:

R2: initially increases with decreasing if is held fixed, asymptotically approaching a constant value as . This result was previously reported by FH16.

R3: R1 and R2 still apply if the conductive heat flux through the top boundary is (partially) replaced by a fixed heat flux designed to mimic smallscale surface convection.
Note that the range of spanned by the simulations presented here is insufficient to provide a robust characterization of the trend found in R1. In particular, the dependence appears roughly linear (Fig. 7), but we cannot rule out other power law exponents or functional relationships. Similar caveats apply to R3. By contrast, R2 has indeed been thoroughly and robustly characterized by FH16.
5.1 RayleighBénard (RB) Convection
To provide context for understanding results R1 and R2 it is instructive to consider the most wellstudied manifestation of thermal convection in the literature, namely RayleighBénard (RB) convection. The physical system consists of a Boussinesq fluid (no density stratification) between two semiinfinite horizontal plates at fixed temperature (Chandrasekhar, 1961; Gastine et al., 2015). The mathematical system can be completely specified by two nondimensional numbers, namely the Rayleigh number and the Prandtl number, . In the Boussinesq limit, , where is the thermal expansion coefficient, is the gravitational acceleration, is the distance between the plates, and is the imposed temperature difference.
The characteristic convective velocity scale, , is typically expressed nondimensionally as the Reynolds number . Whereas and are system parameters, is a diagnostic output of the numerical or laboratory experiment and, for , it is found to scale as
(18) 
If the thermal boundary layer is thinner than the viscous boundary layer () and if it is assumed that there is a local force balance between buoyancy and advection in each convective plume, then the basic mixing length theory suggests and (Siggia, 1994). However, more sophisticated theories and laboratory experiments paint a more complex picture, with multiple parameter regimes determined by the spatial distribution of the thermal and viscous dissipation and the relative widths of the corresponding boundary layers (Ahlers et al., 2009; Gastine et al., 2015). Still, in all of these regimes, and , implying that increases with increasing and decreases with increasing .
This qualitative scaling is straightforward to interpret. For fixed , increasing by, for example, increasing the temperature difference or decreasing the diffusion, will cause the flow to become more turbulent (larger ). Meanwhile, increasing while keeping fixed implies that the flow becomes more viscous (larger ) and laminar (smaller ). However, the dimensional scaling is more subtle. Dimensionally, eq. (18) implies
(19) 
So, for fixed , decreases with decreasing (R1) only if . This is indeed the case for the mixing length estimate quoted above and for most of the parameter regimes described by Ahlers et al. (2009). Meanwhile, becomes independent of at fixed (R2) if . The above estimate of is less than so in this scaling regime, would decrease with decreasing , in contrast with result R2. However, for and , the more comprehensive analysis by Ahlers et al. (2009) suggests values of that are indeed equal to . So the scaling laws quoted in the literature for RB convection are at least qualitatively consistent with our results R1 and R2.
However, a more careful consideration of the underlying physics suggests that the trends we are seeing (relevant for solar convection) are fundamentally different that the trends previously found for RB convection. In RB convection, the decrease of with decreasing for fixed (R1) is accompanied by a decrease in the heat flux through the layer, as quantified by the Nusselt number (Siggia, 1994; Ahlers et al., 2009). Conversely, an increase in is generally accompanied by an increase in the heat flux. By contrast, the heat flux in our simulations, and in the Sun, is fixed by the solar luminosity. We will address this fundamental difference and its implications in §5.2.
In the RB case, the decrease of for decreasing at fixed (R1) can be attributed to what we concisely refer to hereafter as a shift in the power spectrum. In RB convection, as in the simulations presented here, the characteristic width of convective plumes is roughly determined by the thickness of the thermal boundary layer which scales approximately as (Ahlers et al., 2009, FH16). As decreases, the thermal boundary layer becomes thinner and the plumes become narrower, shifting the entire power spectrum toward higher wavenumber. At fixed , this enhances the viscous dissipation. It also enhances the turbulent entrainment of surrounding fluid, suppressing the buoyant driving. The result is a decrease in both (R1) and the heat flux, .
However, if is decreased in pace with (fixed ), then the viscous dissipation can be held in check. Now one might expect to be determined by the potential energy available in the upper boundary layer; see FH16 and §5.4 below. Note that the shift in the power spectrum and in the buoyancy work toward higher wavenumbers as is decreased has been demonstrated by FH16. However, since this is achieved by means of a fixed path through parameter space, this shift is not associated with a significant increase in the viscous dissipation.
5.2 The Luminosity Constraint
In section 5.1 we identified one possible explanation for result R1, namely a shift in the power spectrum toward higher wavenumber that leads to enhanced viscous dissipation for a fixed value of . Here we identify two more possible explanations which refer to concisely hereafter as the Peclet number effect and the raising of the adiabat. Both of these are linked to the requirement that the convection carry a fixed solar luminosity. As such, they have no counterpart in classical RB convection where the convective heat flux varies in general with the dissipation (§5.1).
In our numerical experiments, as in the Sun, a fixed heat flux enters the convection zone from below by means of radiative diffusion. We express this as a radiative heat flux that is determined by the heating function ( is equal to the divergence of ; see eq. 12). This radiative heat flux is the same for each simulation and sets the value of the luminosity that passes through the entire computational domain in equilibrium.
In the mid CZ, the solar luminosity is carried mainly by the convective enthalpy flux, (see Figs. 1 and 6 and eq. 13). Here may be regarded as the temperature deficit of a downward plume relative to its surroundings. This is because the dominant source of thermal fluctuations is the upper thermal boundary layer (as we will argue in the ensuing discussion in this section). If these downward plumes are in pressure equilibrium with their surroundings, then , where is the entropy deficit in the plume. In particular, we can set , where is the mean entropy gradient and is the specific entropy of a typical cool, downflow plume. Both and will in general depend on radius but in what follows we will be concerned with the value of in the mid CZ. Thus, if we equate with and if require that the convective enthalpy flux must transport the solar luminosity through the mid CZ (neglecting the relatively small contribution from the inward kinetic energy flux), then we must have
(20) 
Since is the same for every simulation at all times, must scale inversely with .
Thus, in order to understand results R1, R2, and R3, we must understand the factors that determine . To this end, we consider the curves shown in Fig. 1 (right). Given our fixed entropy upper boundary condition in the C simulations, downward plumes that emanate from the upper boundary layer will be endowed with a specific entropy of relative to the background state. If they retain this specific entropy endowment as they travel downward into the mid CZ, then . Meanwhile, the background entropy in the mid CZ is approximately equal to the entropy contrast across the entire CZ, . This is due both to the efficient mixing of entropy by the convection, which minimizes apart from the upper boundary layer, and to the lower boundary condition of , which can be used to define the base of the CZ
Now the significance of the downflow plumes in determining (and ) becomes clear. A plume that begins in the upper boundary layer and travels downward to the mid CZ will have an entropy deficit that is comparable to . By comparison, fluid flowing up from the lower boundary is nearly isentropic () and most of it overturns before it ever reaches the mid CZ as a consequence of the density stratification. This a characteristic feature of compressible convection (Spruit et al., 1990; Miesch, 2005; Nordlund et al., 2009).
We can generalize this picture by writing , where is a number between zero and one that quantifies how well downflow plumes can hold on to the specific entropy deficit bequeathed to them in the upper boundary layer. For example, if is large, then the thermal content of plumes will diffuse away before they make it to the mid CZ. This would yield a value of close to zero and inefficient convective heat transport regardless of the value of . The relevant nondimensional number here is the Peclet number, which we define here as , where is the width of the thermal boundary layer (equal to the horizontal width of a plume), and is the depth of the CZ. In the limit , plumes will be able to retain their thermal variations with negligible diffusive losses. However, this does not necessarily imply values of close to unity. Other factors also contribute to , most notably the density stratification. As plumes travel downward, they must entrain surrounding fluid in order to maintain their downward momentum. This dilutes their entropy content and introduces a multiplicative contribution to that is inversely proportional to the density contrast between the upper and mid CZ. The value of may also depend on the Reynolds number, which will impact the stability of plumes and may regulate turbulent entrainment.
For our purposes here, we will neglect these additional factors and refer to this phenomena concisely as the Peclet number effect. In short, high values of can promote low velocities by enhancing the thermal content of plumes, as reflected by an increase in the value of . The contribution of the density stratification to is important in general but not relevant here since all simulations have the same density contrast.
If the SGS heat flux in the surface layers is assumed to be diffusive in nature, as in §3, then the Peclet number effect should cease to play a role. This is because the thermal boundary layer in this case scales as , which would make the diffusion time across a plume, independent of . As the diffusion decreases, plumes become smaller, so their diffusive losses are undiminished and their ability to carry heat across the CZ depends only on their speed, . Lower values of would then imply lower values of , which would in turn lead to lower thermal variations . This situation is not consistent with the luminosity constraint (20). However, the Peclet number effect likely does play a role if the width of the boundary layer, , is determined by processes other than conduction, as in §4. We will return to this issue in §5.4.
As noted in §2.2, increases for the C series of simulations as is decreased. This is a consequence of the requirement that conduction must carry the heat flux through the upper boundary, which in turn requires that the upper entropy gradient be (see eq. 16). Since the width of the thermal boundary layer, scales as , this implies . Given our thermal boundary conditions ( at top, at bottom), this effectively raises the specific entropy of the entire CZ, apart from the upper boundary layer.
We refer to this effect as the raising of the adiabat. In short, a decrease in will induce a decrease in by enhancing the thermal deficit of downflow plumes relative to their surroundings, (see eq. 20). It accomplishes this by raising the adiabat of the CZ, enhancing .
5.3 High Prandtl Number Convection
The path in parameter space in which is held fixed while is decreased (result R1) implies increasingly high values of the Prandtl number, . Thermal convection in the regime and has been extensively studied in the literature within the context of mantle convection in terrestrial planets. Notable examples include global simulations of mantle convection that exhibit thin, coherent “superplumes” that span the entire convection zone and dominate the heat transport (Bercovici et al., 1992; Deschamps et al., 2010, 2012; Kronbichler et al., 2012). Similar structures have also been observed in laboratory experiments (Davaille and Limare, 2007). Furthermore, there is some evidence from seismic imaging that such superplumes do in fact exist in the Earth’s mantle (e.g. Romanowicz and Gung, 2002).
It has been suggested by Spruit (1997) and Hanasoge et al. (2016) that thin, coherent plumes roughly analogous to mantle superplumes may also exist in the solar convection zone, and may help resolve the convective conundrum. This is an intriguing possibility that warrants further study. It would be roughly consistent with the scenarios advocated here, in which heat transport is dominated by coherent plumes characterized by high Peclet numbers and thus large thermal contrasts relative to the background (cf. eq. 20).
However, it must be remembered that mantle convection is far different from stellar convection. In particular, the flow is nearly incompressible and highly viscous, characterized by Reynolds numbers much less than unity (in stars ). Furthermore, large variations in viscosity, thermal conductivity, and the thermal expansion coefficient are all thought to play an important role in the formation of terrestrial superplumes (Zhang et al., 1997; Yuen et al., 2007).
5.4 Assessment
We opened this section (§5) with the motivation of understanding result R1; why does decrease with decreasing if is held fixed? In §5.1–5.2 we offered three potential reasons: (1) a shift in the power spectrum (toward higher wavenumbers), (2) the Peclet number effect, and (3) the raising of the adiabat. Only the first applies to the wellstudied problem of RB convection, the other two arising from the constraint that the convection must carry out the solar luminosity.
Which effect is dominant? This is addressed by Figure 8, which shows the distribution of entropy fluctuations in the mid CZ for several cases from section 3 on the constant track through parameter space. This demonstrates an increase in the entropy deficit in cool downflows, , as is decreased. Since (2) the Peclet number effect is expected to be ineffective for these cases with a conductive SGS heat flux (see §5.2), we attribute the decrease in with decreasing (R1) mainly to (3) the raising of the adiabat.
Then what of R2? Why does the raising of the adiabat become ineffective for the fixed path through parameter space? We attribute this to the lower value of , which may destabilize plumes and increase turbulent entrainment. Recall from section 3 that the effective of the flow increases with decreasing on this path, in contrast to the fixed path where the laminarization of plumes (lower ) may promote more efficient heat transport. As plumes lose coherence on the fixed path, the kinetic energy imparted to convection is determined not by the thermal content of individual flow structures, but by the net potential energy available to sustain the convection, which scales as (FH16). As increases (), the thermal boundary layer thins () so that the potential energy available to sustain convection becomes independent of . The result is that becomes independent of dissipation, as was demonstrated by FH16 for .
It may be argued that effects (3) and possibly (1) are artificial in that they depend on our representation of unresolved nearsurface convection as a conductive heat flux proportional to . However, it is clear that our simulations do underestimate and overestimate the thickness of the thermal boundary layer, , relative to the actual Sun. So, we expect that the trends identified here will still be relevant even if the nature of the subgridscale heat transport is nondiffusive.
This expectation is borne out by the results presented in §4 (result R3). Here we partially replaced the conductive SGS heat flux with a fixed heat flux that is independent of both and , in the limit of small . The results were similar to the conductive case (results R1 and R2). However, this cannot be attributed to the raising of the adiabat because in these G cases does not increase with decreasing (Fig. 6). Yet, the thermal fluctuations are still enhanced in the low limit, as shown in Fig. 8. This points to (2) the Peclet number effect as the likely culprit. As in the C series, these effects are suppressed in the constant path through parameter space as the flow becomes more turbulent (higher ) and plumes lose coherence.
We intend to explore these conclusions more thoroughly in the future with higherresolution simulations that achieve low enough values of that the conductive heat flux through the outer surface becomes negligible. In this regime, effects (1) and (3), namely the shift of the power spectrum and the raising of the adiabat, should cease to play a role because they depend on the structure of the boundary layer, which will become essentially independent of . However, the Peclet number effect (2) should become more important in this regime, as noted in §5.2. This is because the width of the plumes should become independent of and should scale as .
It is worth emphasizing explicitly that all three effects should saturate at very low values of , approaching a regime where becomes independent of dissipation ( and ). For (2) the Peclet number effect, this occurs when and thermal diffusion across plumes becomes negligible. For (1) the shifting of the power spectrum and (3) the raising of the adiabat, saturation should occur when approaches the actual depth of the photospheric boundary layer in the Sun. Simulations of surface convection suggest that this is on the order of 100 km, if defined as where the buoyancy driving occurs (Nordlund et al., 2009). However, may be established over a slightly wider region, possibly as much as 1 Mm (e.g. Lord et al., 2014). This suggests that, is only about 12 orders of magnitude thinner than in the highestresolution simulations considered here. Thus, extrapolating the scaling relationship six orders of magnitude or more in would yield boundary layer widths that are unrealistically thin (FH16).
It is also worth noting that in any asymptotic regime in which becomes independent of dissipation, then the mean entropy deficit in plumes, will also become independent of dissipation. This follows from eq. (20).
6 Summary and Implications
In this paper we have argued that, if solar convection can operate in a high Prandtl number regime, then this may help alleviate the convective conundrum (Sec. 1). This has also been argued by Hotta et al. (2015b).
In particular, we have conducted a series of nonrotating, global solar convection simulations with progressively smaller values of the thermal diffusion, , in order to approach the solar parameter regime. If is held fixed as we decrease (corresponding to ), we find that the amplitude of the convective velocity systematically decreases (results R1, R3 in §5). For the cases considered in Section 3 (), we attribute this to an increase in the entropy contrast across the CZ, , which enhances the entropy deficit in cool, downflow plumes. Larger entropy variations in turn implies that the convection can carry the imposed solar luminosity with smaller flow speeds. We expect this behavior to saturate when the true of the solar convection zone is reached, after which should become independent of .
An increase in the effective Peclet number can also help reduce by enhancing thermal fluctuations in the mid CZ. However, we argue that this should only play a role if the modeled SGS heat flux associated with smallscale surface convection is nondiffusive (§5.4). For a diffusive (conductive) heat flux, the width of plumes becomes narrower as is decreased so diffusive losses remain undiminished as the thermal diffusivity is reduced. For a nondiffusive SGS heat flux, we expect that will become independent of when . These mechanisms for regulating are distinct from the wellstudied problem of RayleighBénard convection because they are intimately linked to the requirement that the convection must carry out the solar luminosity (§5.2).
We find that if we decrease as we decrease (fixing ), the results are less favorable from the perspective of the convection conundrum, but still promising. In particular, as is decreased, asymptotes to a constant value that is independent of (result R2 of §5), as found in the more comprehensive study of FH16.
We emphasize again that the simulations reported here are nonrotating. Further work needs to be done to determine whether our results persist for simulations that include rotation or whether other mechanisms for regulating might be present. If they do persist, this raises the question of: Why might solar convection operate in a high regime?
A plausible answer is smallscale magnetism. Hotta et al. (2015b) have recently shown that magnetic fields generated by smallscale dynamo action in highresolution solar convection simulations can simultaneously inhibit shear and reduce turbulent mixing of entropy fluctuations in downflow plumes. In other words, the effective viscosity is enhanced and the effective thermal diffusion is reduced, implying an effective value of that is greater than one. Furthermore, as the resolution is increased, Hotta et al. (2015b) argue that the smallscale magnetic energy should saturate, approaching a value that is close to the equipartition value associated with the kinetic energy of the convection. Thus, one might expect the effective to remain constant as the explicit SGS thermal diffusion (represented by ) is decreased. This provides some justification for our fixed path through parameter space.
The conclusion that smallscale magnetism can act as an effective viscosity is also supported by a number of global convection simulations (Brun et al., 2004; Nelson et al., 2013; Fan and Fang, 2014; Karak et al., 2015). In particular, Fan and Fang (2014) and Karak et al. (2015) demonstrate that the suppression of smallscale shear by magnetism can modify the transition from solar to antisolar differential rotation discussed in Sec. 1, shifting it toward higher Rossby numbers. In other words, the inclusion of magnetism can promote solarlike differential rotation in a simulation that would otherwise be antisolar.
Though promising, these results must still be considered speculative; more work is needed to determine whether or not solar convection effectively operates in the high regime and if so, whether or not this can solve the convection conundrum.
We thank Ben Brown, Sacha Brun, Brad Hindman, Hideyuki Hotta, Bidya Karak, Mark Rast, and Matthias Rempel for many enlightening discussions on the issues covered here and we thank the anonymous referees for helpful comments on the manuscript. This research is supported by NASA grant NNX13AG18G. O’Mara was supported by the Research Experience for Undergraduates (REU) program at University of Colorado’s Laboratory for Atmospheric and Space Physics (LASP) and by a research stipend from Regis University. Featherstone and the development of Rayleigh were supported by the Computational Infrastructure for Geodynamics (CIG), which is supported by the National Science Foundation award NSF094946. This work used computational resources provided through NASA HEC support on the Pleiades supercomputer in conjunction with NASA grant NNX14AC05G. NCAR is sponsored by the National Science Foundation.
Footnotes
 journal: Advances in Space Research
 The base of the CZ may also be defined as the radius where the convective heat flux crosses zero. With this alternative definition, the point where may not precisely coincide with the base of the CZ but this distinction is irrelevant here.
References
 Ahlers, G., Grossmann, S., Lohse, D., 2009. Heat transfer and large scale dynamics in turbulent rayleighbénard convection. Reviews of Modern Physics 81, 503–537.
 Augustson, K.C., Brun, A.S., Miesch, M.S., Toomre, J., 2015. Grand minima and equatorward propagation in a cycling stellar convective dynamo. ApJ 809, 149 (25pp).
 Bercovici, D., Schubert, G., Glatzmaier, G.A., 1992. Threedimensional convection of an infiniteprandtlnumber compressible fluid in a basally heated spherical shell. J. Fluid Mech. 239, 683–719.
 Brun, A.S., Miesch, M.S., Toomre, J., 2004. Globalscale turbulent convection and magnetic dynamo action in the solar envelope. ApJ 614, 1073–1098.
 Chandrasekhar, S., 1961. Hydrodynamic and Hydromagnetic Stability. Oxford University Press, Oxford.
 Chapman, C.J., Proctor, M.R.E., 1980. Nonlinear rayleigh bénard convection between poorly conducting boundaries. J. Fluid Mech. 101, 759–782.
 Charbonneau, P., 2010. Dynamo models of the solar cycle. Living Reviews in Solar Physics 7. Http://www.livingreviews.org/lrsp20103.
 Clune, T.C., Elliott, J.R., Miesch, M.S., Toomre, J., Glatzmaier, G.A., 1999. Computational aspects of a code to study rotating turbulent convection in spherical shells. Parallel Computing 25, 361–380.
 Davaille, A., Limare, A., 2007. Laboratory Studies of Mantle Convection. Treatise on Geophysics, Volume 7: Mantle Dynamics, second edition, chapter 3. Elsevier, Amsterdam. pp. 89–165
 Depassier, M.C., Spiegel, E.A., 1981. The largescale structure of compressible convection. Astronom. J. 86, 496–512.
 Deschamps, F., Tackley, P.J., Nakagawa, T., 2010. Temperature and heat flux scalings for isoviscous thermal convection in spherical geometry. Geophys. J. Int. 182, 137–154.
 Deschamps, F., Yao, C., Tackley, P.J., SanchezValle, C., 2012. High rayleigh number thermal convection in volumetrically heated spherical shells. J. Geophys. Res. 117, E09006.
 Fan, Y., Fang, F., 2014. A simulation of convective dynamo in the solar convective envelope: Maintenance of the solarlike differential rotation and emerging flux. ApJ 789, 35 (13pp).
 Featherstone, N.A., Hindman, B.W., 2016. The spectral amplitude of stellar convection and its scaling in the highrayleighnumber regime. Astrophys. J., submitted.
 Featherstone, N.A., Miesch, M.S., 2015. Meridional circulation in solar and stellar convection zones. ApJ 804, 67 (21pp).
 Gastine, T., Wicht, J., Aurnou, J.M., 2015. Turbulent rayleighbénard convection in spherical shells. J. Fluid Mech. 778, 721–764.
 Gastine, T., Yadav, R.K., Morin, J., Reiners, A., Wicht, J., 2014. From solarlike to antisolar differential rotation in cool stars. MNRAS 438, L76–L80.
 Gilman, P.A., 1983. Dynamically consistent nonlinear dynamos driven by convection in a rotating spherical shell ii. dynamos with cycles and strong feedbacks. ApJSS 53, 243–268.
 Gilman, P.A., Glatzmaier, G.A., 1981. Compressible convection in a rotating spherical shell. i. anelastic equations. ApJSS 45, 335–349.
 Glatzmaier, G.A., 1985. Numerical simulations of stellar convective dynamos. ii. field propogation in the convection zone. ApJ 291, 300–307.
 Gough, D.O., 1969. The anelastic approximation for thermal convection. J. Atmos. Sci. 26, 448–456.
 Greer, B.J., Hindman, B.W., Featherstone, N.A., Toomre, J., 2015. Helioseismic imaging of fast convective flows throughout the nearsurface shear layer. ApJ 803, L17 (5pp).
 Guerrero, G., Smolarkiewicz, P.K., Kosovichev, A.G., Mansour, N.N., 2013. Differential rotation in solarlike stars from global simulations. ApJ 779, 176 (13pp).
 Hanasoge, S.M., Duvall, T., DeRosa, M.L., 2010. Seismic constraints on interior solar convection. ApJL 712, L98–L102.
 Hanasoge, S.M., Duvall, T., Sreenivasan, K.R., 2012. Anamolously weak solar convection. Proc. Nat. Acad. Sci. 109, 11928–11932, doi:10.1073/pnas.1206570109.
 Hanasoge, S.M., Gizon, L., Sreenivasan, K.R., 2016. Seismic sounding of convection in the sun. Ann. Rev. Fluid Mech. 48, 191–217.
 Hotta, H., Rempel, M., Yokoyama, T., 2014. Highresolution calculations of the solar global convection with the reduced speed of sound technique. i. the structure of the convection and the magnetic field without the rotation. ApJ 786, 24 (18pp).
 Hotta, H., Rempel, M., Yokoyama, T., 2015a. Highresolution calculations of the solar global convection with the reduced speed of sound technique. ii. near surface shear layer with the rotation. ApJ 798, 51 (15pp).
 Hotta, H., Rempel, M., Yokoyama, T., 2015b. Efficient smallscale dynamo action in the solar convection zone. ApJ 803, 42 (14pp).
 Jones, C.A., Boronski, P., Brun, A.S., Glatzmaier, G.A., Gastine, T., Miesch, M.S., Wicht, J., 2011. Anelastic convectiondriven dynamo benchmarks. Icarus 216, 120–135.
 Käpylä, P.J., Käpylä, M.J., Brandenburg, A., 2014. Confirmation of bistable stellar differential rotation profiles. A&A 570, A43 (10pp).
 Karak, B.B., Käpylä, P.J., Käpylä, M.J., Brandenburg, A., Olspert, N., Pelt, J., 2015. Magnetically controlled stellar differential rotation near the transition from solar to antisolar profiles. A&A 576, A26 (17pp).
 Kronbichler, M., Heister, T., Bangerth, W., 2012. High accuracy mantle convection simulation through modern numerical methods. Geophys. J. Int. 191, 12–29.
 Lord, J., Cameron, R., Rast, M., Rempel, M., Roudier, T., 2014. The role of subsurface flows in solar surface convection: Modeling the spectrum of supergranular and larger scale flows. ApJ 793, 24 (11pp).
 Miesch, M.S., 2005. Largescale dynamics of the convection zone and tachocline. Living Reviews in Solar Physics 2. Http://www.livingreviews.org/lrsp20051.
 Miesch, M.S., Featherstone, N.A., Rempel, M., Trampedach, R., 2012. On the amplitude of convective velocities in the deep solar interior. ApJ 757, 128 (14pp).
 Nelson, N.J., Brown, B.P., Brun, A.S., Miesch, M.S., Toomre, J., 2013. Magnetic wreathes and cycles in convective dynamos. ApJ 762, 73 (20pp).
 Nordlund, A., Stein, R.F., Asplund, M., 2009. Solar surface convection. Living Reviews in Solar Physics 6. Http://www.livingreviews.org/lrsp20092.
 Racine, E., Charbonneau, P., Ghizaru, M., Bouchat, A., Smolarkiewicz, P.K., 2011. On the mode of dynamo action in a global largeeddy simulation of solar convection. ApJ 735, 46 (22pp).
 Romanowicz, B., Gung, Y., 2002. Superplumes from the coremantle boundary to the lithosphere: Implications for heat flux. Science 296, 513–516.
 Siggia, E.D., 1994. High rayleigh number convection. Annu. Rev. Fluid Mech. 26, 137–168.
 Spruit, H.C., 1997. Convection in stellar envelopes: A changing paradigm. Mem. Soc. Astr. It. 68, 397–414.
 Spruit, H.C., Nordlund, A., Title, A.M., 1990. Solar convection. Annu. Rev. Astron. Astrophys. 28, 263–301.
 Yuen, D.A., Monnereau, M., Hansen, U., Kameyama, M., Matyska, C., 2007. Dynamics of superplumes in the lower mantle, in: Yuen, D.A., Maruyama, S., Karato, S.I., Windley, B.F. (Eds.), Superplumes: Beyond Plate Tectonics, Springer, Dordrecht. pp. 239–267.
 Zhang, J., Childress, S., Libchaber, A., 1997. Nonboussinesq effect: Thermal convection with broken symmetry. Phys. Fluids 9, 1034–1042.