Prandtl-number Effects in High-Rayleigh-number Spherical Convection

The Prandtl number, $\Pr$ = ν/κ, where ν is the kinematic viscosity and κ is the thermal diffusivity, is known to be small in the dynamo regions of planets and stars (Ossendrijver 2003; Roberts 2007). It is well known from linear theory that the value of the Prandtl number both influences the structure and amplitude of convective motions and controls the critical Rayleigh number required for the onset of rotating convection (e.g., Chandrasekhar 1961). Asymptotic approximations to the governing fluid equations can be carried out based on the size of the Prandtl number. Such approximations yield some insight into the convective dynamics arising under different Prandtl-number regimes. In the large Prandtl-number limit, the influence of inertia is weak and can be neglected; this limit is routinely exploited in studying the convection of planetary mantles where the Prandtl number can be of order 10²⁰ (Schubert et al. 2001). Spiegel (1962) developed an approximate set of equations valid in the limit $\Pr \to 0$ that showed that inertia plays a leading-order role in the convective dynamics (see also Thual 1992). Both of these approximate models have been applied only to nonrotating and incompressible Boussinesq systems to date.

Much less is known about the role of the Prandtl number in compressible convection. Earlier studies of compressible convection have primarily used a Prandtl number of order unity (e.g., Gilman 1977; Gilman & Glatzmaier 1981; Goudard & Dormy 2008; Christensen 2011; Schrinner et al. 2012; Soderlund et al. 2012; Gastine et al. 2015, 2016; Wicht & Meduri 2016). While there are also many studies that make use of nonunity Prandtl numbers (e.g., Brown et al. 2011; Käpylä et al. 2013; Jones 2014; Nelson et al. 2014; Augustson et al. 2015, 2016; Duarte et al. 2016; Brun et al. 2017), no parameter studies that vary the Prandtl number systematically have been carried out. One exception to this trend is the work of O'Mara et al. (2016), who explored the characteristics of high Prandtl number compressible convection. These authors found that high Prandtl number convection tended to possess lower characteristic flow speeds with respect to unity Prandtl-number convection, owing to the enhanced entropy content of its downflow plumes. Finally, we note that recent work using a small Prandtl number with rapid rotation has found that the anelastic approximation can yield spurious behavior (Calkins et al. 2015a) that does not appear in nonrotating anelastic convection (Calkins et al. 2015b). These results raise serious questions regarding the applicability of the anelastic approximation within rotating stellar interiors where it remains to be seen if the convective flows are well-approximated by a Prandtl number of unity.

Systematic parameter space studies of convective dynamics in stellar interiors have, so far, focused largely on the role of buoyancy driving and rotation. Featherstone & Hindman (2016b) investigated the response of the convection to varying Rayleigh numbers and varying degrees of density stratification. Those simulations were nonrotating, hydrodynamic, $\Pr$ = 1, and demonstrated a clear scaling relationship between kinetic energy and Rayleigh number. Those results also suggest that a naive interpretation of model results (by ascribing solar values to all of the problem parameters but the diffusion coefficients) will naturally overestimate the low-wavenumber power in the convective power spectrum. The influence of rotation was investigated using a similar methology by Featherstone & Hindman (2016a) who identified a complementary scaling law relating convective-cell size and rotational influence. When rotation is present, and diffusive effects are negligible, the typical spatial size of convective cells is determined primarily through the Rossby number, which expresses the ratio of the rotation period to a characteristic convective timescale. Their work was also restricted to $\Pr$ = 1.

Through this paper, we extend these studies and examine the effects of Prandtl-number variation on the convective dynamics. We present a series of numerical simulations designed to examine how the structure and amplitude of the convective flow within a stellar interior depends on the Prandtl number and the convective forcing. We vary the Prandtl number and the convective forcing in a systematic way, covering both low and high Prandtl numbers. Effects due to rotation and magnetism are not included. We will show that the convection develops smaller-scale structures as the convective forcing is increased and the Prandtl number is decreased, corresponding to an increase in high-wavenumber power. As the high-wavenumber power increases, the low-wavenumber power decreases and this trend occurs for all of the Prandtl numbers studied. We also show that the Prandtl number has an important influence on the boundary-layer thickness.

We begin our examination of the convective energetics by looking at the integrated kinetic energy KE, defined as

$$\begin{eqnarray}&&\mathrm{KE}=\displaystyle \frac{1}{2}\int \bar{\rho }(r){| {\boldsymbol{u}}(r,\theta ,\phi )| }^{2}\ {d}^{3}x,\end{eqnarray} \tag{ 11 }$$

where the integration is computed over the entire domain and then time averaged.

Figure 1(a) shows the dimensional kinetic energy versus ${\mathrm{Ra}}_{F}$. Different symbols indicate different values of $\Pr$. As the Prandtl number is lowered, the kinetic energy increases for any given ${\mathrm{Ra}}_{F}$. The kinetic energy for those runs with ${\mathrm{Ra}}_{F}$ ≳ 2 × 10⁴ appears to have reached a steady value that is independent of ${\mathrm{Ra}}_{F}$. The level at which the kinetic energy saturates is dependent on the Prandtl number. The saturation of the kinetic energy as the Rayleigh number is changed was also found in Featherstone & Hindman (2016b), although their study was restricted to a Prandtl number of unity.

Zoom In Zoom Out Reset image size

As $\Pr$ is increased, ν becomes larger leading to enhanced viscous dissipation and a smaller Reynolds number. The Reynolds number measures the relative ratio of inertial forces to viscous forces and is given by

$$\begin{eqnarray}&&\mathrm{Re}=\displaystyle \frac{\sqrt{\widetilde{{| {\boldsymbol{u}}| }^{2}}}H}{\nu },\end{eqnarray} \tag{ 12 }$$

where H and the tilde retain the same meaning as before, representing the shell depth and a volume average, respectively. Larger Prandtl numbers will produce smaller Reynolds numbers for a given ${\mathrm{Ra}}_{F}$.

Figure 1(a) indicates that below some ${\mathrm{Ra}}_{F}$ cutoff, diffusion plays a leading-order role in the force balance. Beyond ${\mathrm{Ra}}_{F}$ ∼ 2 × 10⁴, which we denote the high-${\mathrm{Ra}}_{F}$ regime,⁶ this is no longer the case; diffusion no longer plays a leading-order role in the global force balance and the kinetic energy remains constant as ${\mathrm{Ra}}_{F}$ is increased. We note that the cutoff for the high-${\mathrm{Ra}}_{F}$ region is based on the $\Pr$ = 1 results of Featherstone & Hindman (2016b). Importantly, the cutoff is a decreasing function of the Prandtl number; lower Prandtl numbers will have a lower high-${\mathrm{Ra}}_{F}$ cutoff. To simplify the analysis, only a single cutoff is used.

The kinetic energy can also be discussed from a nondimensional point of view. To do this, we choose a nondimensional measure of the kinetic energy $\widehat{\mathrm{KE}}$ as

$$\begin{eqnarray}&&\widehat{\mathrm{KE}}\equiv \displaystyle \frac{{H}^{2}}{{\nu }^{2}M}\mathrm{KE}.\end{eqnarray} \tag{ 13 }$$

The nondimensionalization has been carried out using the mass M contained within the spherical shell, the shell depth, and the viscous diffusion timescale. Under such a nondimensionalization, the kinetic energy can be related to the Reynolds number as

$$\begin{eqnarray}&&\widehat{\mathrm{KE}}\sim {\mathrm{Re}}^{2}.\end{eqnarray} \tag{ 14 }$$

Therefore we do not plot the nondimensional kinetic energy separately, but do list it for each run in Table 1.

If ${\mathrm{Ra}}_{F}$ is large enough to be in the high-${\mathrm{Ra}}_{F}$ regime, then the kinetic energy is largely insensitive to the level of diffusion. This implies that the velocity is not strongly dependent on the level of diffusion and that the Reynolds number should scale with the viscosity as $\mathrm{Re}$ ∼ ν⁻¹. Given our definition of ${\mathrm{Ra}}_{F}$ and $\Pr$, namely that ${\mathrm{Ra}}_{F}$ ∼ ν⁻¹κ⁻² and $\Pr$ ∼ νκ⁻¹, the Reynolds number scaling becomes $\mathrm{Re}\,\sim \,{\mathrm{Ra}}_{F}^{1/3}{\Pr }^{-2/3}$. In Figure 1(b), we plot the Reynolds number versus ${\mathrm{Ra}}_{F}$ $\Pr$⁻², where the different symbols are the same as in panel (a). Each data point is time averaged. A least-squares fit to the data in the high-${\mathrm{Ra}}_{F}$ region yields a scaling law of

$$\begin{eqnarray}&&\mathrm{Re}\,\propto \,{({\mathrm{Ra}}_{F}{\Pr }^{-2})}^{0.373\pm 0.008}.\end{eqnarray} \tag{ 15 }$$

These results indicate that there are two distinct parameter regimes: one in which diffusion is an important factor in the global force balance (the low-${\mathrm{Ra}}_{F}$ region), and one where diffusion no longer plays an appreciable role in the interior, bulk global force balance (the high-${\mathrm{Ra}}_{F}$ region).

We find that the kinetic energy may reach a ${\mathrm{Ra}}_{F}$-independent regime, but the flow's morphology is still affected by the level of diffusion. This can be seen in the relative spectral distribution of velocity between the high- and low-${\mathrm{Ra}}_{F}$ systems. Figure 2 shows the velocity power spectra for all of the runs with Prandtl numbers of $\Pr$ = 0.25, $\Pr$ = 1.0, and $\Pr$ = 4.0. Each spectrum has been normalized such that it has unit integrated power. The rows correspond to the different Prandtl numbers. The first column shows the spectra taken at the lower convection zone or r/r_o ≈ 0.775. The second column shows the spectra near the upper boundary, or r/r_o ≈ 0.985, in the thermal boundary layer. In each panel, all of the spectra with the given Prandtl number are plotted. Each spectrum is colored by ${\mathrm{Ra}}_{F}$ $\Pr$, with high values taking on red tones and low values displaying blue tones. Each spectrum is a time average over several tens of overturnings.

Zoom In Zoom Out Reset image size

Across all of the cases, as ${\mathrm{Ra}}_{F}$ increases from low values, the point-wise velocity power increases at nearly all spherical harmonic degrees. The ℓ-value associated with the peak for each spectra increases with increasing ${\mathrm{Ra}}_{F}$. This indicates that smaller-scale structures become more apparent with increasing ${\mathrm{Ra}}_{F}$. At sufficiently high Rayleigh number, power in the high-ℓ portion continues to increase, but the low-ℓ portion starts to decrease. This suggests a break down of large-scale coherent structures. The highest Rayleigh number runs have lower power in the large scales compared to the low-${\mathrm{Ra}}_{F}$ counterparts. The low-wavenumber region occurs in the approximate range of ℓ ≲ 10 for most of the simulations in this study. This trend occurs for all values of $\Pr$ that were studied. If viscosity played a significant role in the asymptotic regime, one might expect there to be some variation in the spectra when the Prandtl number is varied. We do not observe large variations between the spectra for the Prandtl numbers within our range indicating that viscosity only plays a minor role in the asymptotic regime.

The third column in Figure 2 plots shell slices of the radial velocity taken near the upper boundary, the same location in radius as the second column. Each slice is a single snapshot in time. The color scale is the same for all three shell slices with red tones indicating positive outward flows and blue tones indicating negative inward flows. The lower Prandtl number run (the top panel) displays larger velocities and more small-scale structures compared to the larger Prandtl cases. The Reynolds numbers of all three slices cover a large range. The top row $\Pr$ = 0.25 run has $\mathrm{Re}$ = 261.2, the middle slice with $\Pr$ = 1.0 has $\mathrm{Re}$ = 56.3, and the bottom slice with $\Pr$ = 4.0 has $\mathrm{Re}$ = 11.7. The large range in Reynolds numbers indicates that the inertial subrange for each simulation is different.

At a sufficiently high Rayleigh number, the integrated KE becomes independent of the level of diffusion (both thermal and momentum diffusion). As ν is decreased and ${\mathrm{Ra}}_{F}$ is increased, the flow becomes more turbulent with smaller-scale structures. To leading-order, once the high-${\mathrm{Ra}}_{F}$ regime is reached, the total integrated dimensional KE is constant as the diffusion is further reduced. This fact is largely independent of $\Pr$ for the Prandtl numbers that were within our range, however, a weak $\Pr$ dependence remains because the inertial subrange is extended as the viscosity is reduced individually. These trends are similar to those found in the simulations of Featherstone & Hindman (2016b), which used a Prandtl number of unity.

The energy transport across the layer can be characterized by four radial energy fluxes: the enthalpy flux F_e, the kinetic energy flux F_KE, the conductive flux F_c, and the viscous flux F_ν, which we define as

$$\begin{eqnarray}&&{F}_{e}=\bar{\rho }{c}_{p}{u}_{r}T,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{KE}}=\tfrac{1}{2}\bar{\rho }{u}_{r}{| {\boldsymbol{u}}| }^{2},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{F}_{c}=\kappa \bar{\rho }\bar{T}\displaystyle \frac{\partial S}{\partial r},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{F}_{\nu }=-({\boldsymbol{u}}\,{\boldsymbol{\cdot }}\,{ \mathcal D })\,{\boldsymbol{\cdot }}\,\hat{r},\end{eqnarray} \tag{ 19 }$$

respectively. Note that the conductive flux F_c is associated with the diffusion of entropy perturbations, and it should not be confused with radiative diffusion arising from the reference state temperature gradient; that effect is represented by Q in our models. Averages are taken of these fluxes over several diffusion times, indicated using brackets. We consider the contribution of conduction by looking at the fractional convective flux ${f}_{\mathrm{conv}}$ defined as

$$\begin{eqnarray}&&{f}_{\mathrm{conv}}\equiv \displaystyle \frac{\int \langle {F}_{e}+{F}_{\nu }+{F}_{\mathrm{KE}}\rangle {d}^{3}x}{\int \langle {F}_{c}+{F}_{e}+{F}_{\nu }+{F}_{\mathrm{KE}}\rangle {d}^{3}x}.\end{eqnarray} \tag{ 20 }$$

Upon rearranging this quantity, we can write

$$\begin{eqnarray}&&\displaystyle \frac{1}{1-{f}_{\mathrm{conv}}}=\displaystyle \frac{\int \langle {F}_{c}+{F}_{e}+{F}_{\nu }+{F}_{\mathrm{KE}}\rangle {d}^{3}x}{\int \langle {F}_{c}\rangle {d}^{3}x}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&=\,1+\displaystyle \frac{\int \langle {F}_{e}+{F}_{\nu }+{F}_{\mathrm{KE}}\rangle {d}^{3}x}{\int \langle {F}_{c}\rangle {d}^{3}x}.\end{eqnarray} \tag{ 22 }$$

This resembles the traditional Nusselt number, but differs in two important ways. First, the traditional Nusselt number does not include the viscous flux or the kinetic energy flux, both of which we include. Second, the conductive flux that appears in our definition is the established conductive flux, not the conductive flux in the absence of convection. Values of $1/(1-{f}_{\mathrm{conv}})$ that are of order unity translate to a lack of convective heat transport. Large values indicate that convection plays a dominant role over thermal conduction in transporting energy through the shell. Figure 3(a) plots $1/(1-{f}_{\mathrm{conv}})$ as a function of ${\mathrm{Ra}}_{F}$ $\Pr$. Viscosity is not expected to play a large role in the heat transport across the shell, which is why we plot against ${\mathrm{Ra}}_{F}$ $\Pr$ and not simply ${\mathrm{Ra}}_{F}$. The plot shows a clear trend that can be fit using least-squares to obtain

$$\begin{eqnarray}&&\displaystyle \frac{1}{1-{f}_{\mathrm{conv}}}\propto {({\mathrm{Ra}}_{F}\Pr )}^{0.275\pm 0.002}\end{eqnarray} \tag{ 23 }$$

with a scaling exponent that is approximately 2/7 (∼0.286). Other studies have found that the Nusselt number scales with the Rayleigh number to the 2/9 power (∼0.222; e.g., Gastine et al. 2015), when a flux based Rayleigh number is used. Our results have a steeper exponent because we include the viscous flux and the kinetic energy flux, which act to increase $1/(1-{f}_{\mathrm{conv}})$ at any given Rayleigh number.

Zoom In Zoom Out Reset image size

We can further characterize the role of conduction by determining how ${\mathrm{Ra}}_{F}$ and $\Pr$ control the thermal boundary-layer thickness. We define the thermal boundary-layer thickness in terms of the time-averaged mean entropy, such that

$$\begin{eqnarray}&&{\delta }_{\mathrm{BL}}=\int \displaystyle \frac{\sup \langle S(r)\rangle -\langle S(r)\rangle }{\sup \langle S(r)\rangle }\ {dr},\end{eqnarray} \tag{ 24 }$$

where the brackets indicate a time-average as before and sup indicates the supremum. Figure 4 shows three examples of the time-averaged mean entropy profiles and the associated boundary-layer location as calculated using Equation (24). We plot the variation of δ_BL with ${\mathrm{Ra}}_{F}$ $\Pr$ for different Prandtl numbers in Figure 3(b). We plot against ${\mathrm{Ra}}_{F}$ $\Pr$ because we do not expect the thermal boundary to depend on the viscosity and the quantity ${\mathrm{Ra}}_{F}$ $\Pr$ is independent of the viscosity and scales as ${\mathrm{Ra}}_{F}$ $\Pr$ ∼ κ⁻³. A least-squares fit gives the scaling law

$$\begin{eqnarray}&&{\delta }_{\mathrm{BL}}\propto {({\mathrm{Ra}}_{F}\Pr )}^{-0.164\pm 0.003}.\end{eqnarray} \tag{ 25 }$$

The scaling exponent is very close to −1/6 (∼ −0.166), indicating that the boundary-layer width scales purely as the thermal diffusivity, i.e., ${\delta }_{\mathrm{BL}}\propto \sqrt{\kappa }$. There is no strong dependence on viscosity within the range of $\Pr$ studied here.

Zoom In Zoom Out Reset image size

We note that the highest Rayleigh number run had a boundary-layer width of about 3% of the shell depth (∼5 Mm). Our thermal boundary layer lacks the radiative processes at work in a star like the Sun, but its physical extent is confined to a similarly small region of the convective domain.

Figure 1 showed that a free-fall regime that is independent of both viscosity and thermal diffusion could be obtained for different Prandtl numbers. We find that in order to have a bulk kinetic energy that is independent of both viscosity and thermal diffusion, the thermal boundary-layer thickness must depend on both the Rayleigh number and the Prandtl number with the scaling law given in Equation (25).

For a fixed Prandtl number, this scaling relation suggests that the boundary-layer thickness will continually decrease as the Rayleigh number is increased. In global simulations that employ explicit diffusivities, the Prandtl number is therefore critical in maintaining a boundary layer that is confined to a small region of the convective domain without becoming vanishingly small as the Rayleigh number is increased.

Certainly, simulations with a conductive boundary layer cannot match every feature of the Sun's thermal boundary layer, for the simple reason that the Sun's boundary layer is regulated by radiative transfer instead of thermal conduction. In a global simulation without radiative transfer, one hopes that the microphysics of the cooling layer can be ignored and only the gross properties of the boundary layer (thickness and entropy contrast) are important. When the boundary layer is conductive, the entropy contrast and thickness are inherently linked. Thus, the best one could hope to do is achieve a physically realistic thickness; this requires that the product $H{({\mathrm{Ra}}_{F}\Pr )}^{-1/6}$ take on the desired thickness, i.e., as the Rayleigh number is increased, the Prandtl number must be decreased. This realistic thickness will probably not, however, coincide with a convective power spectrum that possesses a realistic inertial range.