The Effects of Radially Varying Diffusivities on Stellar Convection Zone Dynamics

Convection is ubiquitous in stellar and planetary interiors where it likely plays an integral role in the generation of magnetic fields. As the interiors of these objects remain hidden from direct observation, numerical models of convection are an important tool in the study of astrophysical dynamos. In such models, unrealistic large values of the viscous ($ν$) and thermal ($κ$) diffusivity are routinely used as an ad-hoc representation of the effects of subgrid scale turbulence which is otherwise too small-scale to resolve numerically. However, the functional forms of these diffusion coefficients can vary greatly between studies, complicating efforts to compare between results and against observations. We explore this issue by considering a series of non-rotating, non-magnetic, solar-like convection models with varying radial functions for the diffusivities and differing boundary conditions. We find that the bulk kinetic energy scales similarly regardless of the diffusivity parameterization. This scaling is consistent with a free-fall scaling, wherein viscosity plays a subdominant role in the force balance. We do not, however, observe such diffusion-free behavior in the convective heat transport. Our results also indicate that the functional form adopted for the diffusion coefficients can impact the distribution of turbulence within the convective shell. These results suggest that some care should be taken when comparing solar convection models directly against helioseismic observations.

💡 Research Summary

This paper investigates how the radial dependence of viscous (ν) and thermal (κ) diffusivities influences the dynamics of a solar‑like convection zone in three‑dimensional, non‑rotating, non‑magnetic simulations. Using the pseudo‑spectral Rayleigh code, the authors model a spherical shell that spans the inner four density scale heights of the Sun, employing the anelastic approximation to capture subsonic, low‑Mach‑number flows. The background stratification follows a polytropic reference state, and the simulations are initialized with a small random entropy perturbation and evolved until a statistically steady state is reached.

The key methodological innovation is the systematic variation of the diffusivity profiles according to ν, κ ∝ ρ^α, where the exponent α takes values –1, –0.5, 0, and 0.5 for internally heated models, while fixed‑entropy and fixed‑flux cases use α = 0. All runs share a Prandtl number of unity, and the primary control parameter is a flux‑based Rayleigh number Ra_F that measures the ratio of buoyancy driving to diffusive damping. A subset of simulations also uses an entropy‑contrast Rayleigh number Ra_S when the entropy difference across the shell is prescribed. Boundary conditions are stress‑free and impenetrable for velocity; entropy or flux conditions are applied at the top and bottom, and an internal heating term Q(r) mimics solar radiative heating.

Across 33 distinct models, the authors examine three principal diagnostics: (1) the total kinetic energy (KE) of the convective flow, (2) the efficiency of heat transport (often expressed as a Nusselt number), and (3) the spectral distribution of kinetic energy as a function of spherical harmonic degree ℓ at several radii.

The first major result is that KE scales with Ra_F in essentially the same way for all α values. At sufficiently high Ra_F the kinetic energy follows a free‑fall scaling (approximately KE ∝ Ra_F^{2/3}), indicating that viscous forces are subdominant in the global force balance. Models with negative α (higher diffusivity near the bottom) reach the asymptotic KE level slightly earlier than those with α = 0 or positive α, but the final KE values converge. This confirms that, for the bulk dynamics, the precise radial profile of ν and κ does not alter the overall vigor of convection.

In contrast, heat transport shows a strong dependence on α. When α > 0, viscosity and thermal diffusivity decrease toward the outer boundary, allowing the upper layers to transport heat efficiently, while the lower layers become more diffusive and thus less effective at carrying flux. For α < 0 the opposite occurs: the lower part of the shell experiences higher viscosity, suppressing convective heat transport there, and the upper part suffers from excessive thermal diffusion, leading to larger conductive losses. Consequently, the Nusselt number and related measures of convective efficiency vary significantly with the chosen diffusivity profile, demonstrating that “diffusion‑free” convection is not realized in the thermal balance even when kinetic energy appears insensitive to diffusion.

Spectral analysis reveals that the radial distribution of turbulent power shifts with α. Positive α concentrates strong, small‑scale motions (high ℓ) in the upper convection zone (r ≈ 0.98 r_o), whereas negative α pushes comparable small‑scale power toward the lower zone (r ≈ 0.75 r_o). The radial velocity snapshots and time‑averaged power spectra illustrate this migration clearly. Thus, the functional form of ν(r) and κ(r) directly reshapes where in the shell turbulence is most intense and at what spatial scales it operates.

The authors also compare internally heated runs with fixed‑entropy and fixed‑flux configurations. Fixed‑entropy cases, governed by Ra_S, display slightly lower KE and heat‑transport efficiencies, while fixed‑flux runs resemble the internally heated models most closely because the imposed flux directly controls the energy budget.

From a broader perspective, the study highlights a subtle but important source of model‑to‑model variability in solar convection simulations. Many existing global convection models adopt large‑eddy diffusivities that are either uniform or follow ad‑hoc radial scalings. This work shows that such choices can bias the depth‑dependent turbulence intensity and the convective heat flux, potentially leading to mismatches with helioseismic inferences of interior flow speeds and temperature gradients. The authors therefore recommend that future comparisons between simulations and observations explicitly state the diffusivity profile used, and, where possible, adopt physically motivated radial dependencies rather than purely numerical convenience.

In conclusion, while the bulk kinetic energy of solar‑like convection appears robust to the details of radial diffusivity, the distribution of turbulence and the efficiency of heat transport are not. Careful treatment of ν(r) and κ(r) is essential for realistic modeling of stellar convection zones and for meaningful interpretation of helioseismic data. The paper suggests that incorporating more realistic subgrid‑scale models—perhaps dynamic Smagorinsky or other LES approaches—could reduce the sensitivity to arbitrary diffusivity choices and bring simulations into closer alignment with observations.