Turbulent Prandtl number from isotropically forced turbulence

Turbulent motions enhance the diffusion of large-scale flows and temperature gradients. Such diffusion is often parameterized by coefficients of turbulent viscosity ($\nu_{\rm t}$) and turbulent thermal diffusivity ($\chi_{\rm t}$) that are analogous to their microscopic counterparts. We compute the turbulent diffusion coefficients by imposing large-scale velocity and temperature gradients on a turbulent flow and measuring the response of the system. We also confirm our results using experiments where the imposed gradients are allowed to decay. To achieve this, we use weakly compressible three-dimensional hydrodynamic simulations of isotropically forced homogeneous turbulence. We find that the turbulent viscosity and thermal diffusion, as well as their ratio the turbulent Prandtl number, ${\rm Pr}_{\rm t} = \nu_{\rm t}/\chi_{\rm t}$, approach asymptotic values at sufficiently high Reynolds and Pecl\'et numbers. We also do not find a significant dependence of ${\rm Pr}_{\rm t}$ on the microscopic Prandtl number ${\rm Pr} = \nu/\chi$. These findings are in stark contrast to results from the $k-\epsilon$ model which suggests that ${\rm Pr}_{\rm t}$ increases monotonically with decreasing ${\rm Pr}$. The current results are relevant for the ongoing debate of, for example, the nature of the turbulent flows in the very low ${\rm Pr}$ regimes of stellar convection zones.


Introduction
The fluids in stellar convection zones are generally characterized by a low microscopic Prandtl number, Pr = / , where is the kinematic viscosity and is the thermal diffusivity (e.g. Ossendrĳver 2003;Augustson et al. 2019). Typical values in the bulk of the solar convection zone, for example, range between 10 −6 and 10 −3 (Schumacher & Sreenivasan 2020). Recently, several studies have explored the possibility that solar convection operates at a high effective Prandtl number regime, meaning that the turbulent Prandtl number Pr t exceeds unity (e.g. O'Mara et al. 2016;Bekki et al. 2017;Karak et al. 2018), as a possible solution to the too high velocity amplitudes in simulations in comparison to the Sun (e.g. Hanasoge et al. 2012;Schumacher & Sreenivasan 2020). However, few attempts have been made to actually measure the turbulent Prandtl number from simulations. A notable exception is the study of Pandey et al. (2021) who reported that the turbulent Prandtl number decreases steeply as a function of the molecular Prandtl number such that Pr t ∝ Pr −1 in simulations of standard Boussinesq and variable heat conductivity Boussinesq convection.
Here we set out to measure turbulent viscosity and thermal diffusivity from a simpler system of isotropically forced homogeneous turbulence. This is done by imposing largescale gradients of velocity and temperature (equivalently specific entropy) and measuring the response of the system. The turbulent diffusion coefficients are computed from the Boussinesq ansatz and an analogous expression for the enthalpy flux. This method provides a direct measurement of the diffusion coefficients without the need to resort to turbulent closures. Similar methods were used recently to measure the turbulent magnetic Prandtl number (Käpylä et al. 2020). We compare our results with those from the widely used − model which was also used by Pandey et al. (2021). We show that the direct results and those from the − model are systematically different and that the latter yields misleading results.

The model
We model isotropically forced, non-isothermal, turbulence in a fully periodic cube of volume (2 ) 3 . We solve the equations of fully compressible hydrodynamics where / = / + · ∇ is the advective derivative, is the density, is the velocity, is the pressure, is the kinematic viscosity, S is the traceless rate-of-strain tensor with is the external forcing, is a relaxation timescale, and 0 is the target mean velocity profile. Furthermore, is the temperature, is the specific entropy, F rad is the radiative flux, C is a cooling term, and 0 is the target mean specific entropy profile. Radiation is modeled via the diffusion approximation, with the radiative flux given by where is the specific heat in constant pressure and is the thermal diffusivity. The ideal gas equation of state = ( − ) = R is assumed, where R is the gas constant, and is the specific heat capacity at constant volume. In the presence of an imposed large-scale flow, viscous dissipation of kinetic energy acts as a source for thermal energy and leads to a linear increase of the temperature. Additional volumetric cooling is applied to counter this with where 0 is the volume-averaged initial temperature and cool is a cooling timescale. We use = cool = ( 0 1 ) −1 , where 0 is the initial uniform value of the sound speed and 1 is the wavenumber corresponding to the box scale.
The external forcing is given by (see Brandenburg 2001) where is the position vector and ( ) = 0 s ( ( ) s / ) 1/2 is a normalization factor where 0 is the forcing amplitude, = | |, is the length of the time step, and − < ( ) < is a random delta-correlated phase. The vector describes non-helical transverse waves, and is given by whereˆis an arbitrary unit vector, and where the wavenumber is randomly chosen. The target profiles of mean velocity and specific entropy are given by (2.10) In addition to the physical diffusion, the advective terms in (2.1) to (2.3) are implemented in terms of fifth-order upwinding derivatives with sixth-order hyperdiffusive corrections and flow-dependent diffusion coefficients; see Appendix B of Dobler et al. (2006). The P C (Pencil Code Collaboration et al. 2021) †, which uses high-order finite differences for spatial and temporal discretisation, was used to produce the numerical simulations.

Units, system parameters, and diagnostics
The equations are non-dimensionalized by choosing the units (2.11) where 0 is the initial uniform density and s0 = √︁ R 0 is the sound speed corresponding to the initial temperature 0 . The level of velocity fluctuations is determined by the forcing amplitude 0 along with the kinematic viscosity. A key system parameter is the ratio of kinematic viscosity and thermal diffusion or the Prandtl number which is varied between 0.01 and 10 in the present study. The Reynolds and Péclet numbers quantify the level of turbulence of the flows: where rms = √︁ ( − 0 ) 2 is the volume-averaged fluctuating rms-velocity and f is the average forcing wavenumber characterizing the energy injection scale. The latter is chosen from a uniformly distributed narrow range in the vicinity of 5 1 . The imposed gradients of large-scale flow and entropy are quantified by where Ma is the Mach number of the mean flow and = = 1 . The Mach number of the turbulent flow is given by where s is the volume-averaged speed of sound. Mean values are taken to be horizontal averages denoted by overbars, that is ( 2.16) Often an additional time average over the statistically steady part of the simulation is taken. Volume averages are denoted by angle brackets . apart from the rms-values which are always assumed to be volume-averaged unless otherwise stated. Errors were estimated by dividing the time series in three parts and averaging over each subinterval. The greatest deviation from the average computed over the whole time series was taken as the error estimate.

Results
The simulations discussed in the present study are listed in table 1.

Turbulent viscosity and heat diffusion from imposed flow and entropy methods
We measure the turbulent viscosity and thermal diffusivity in two ways. First, we impose sinusoidal large-scale profiles of velocity (2.9) or entropy (2.10). The response of the system are non-zero Reynolds stress and vertical enthalpy flux profiles that are parameterized with gradient diffusion terms (e.g. Rüdiger 1989) where primes denote fluctuations from the mean, e.g., = − . The Mach number in the current simulations is always less than 0.1. Therefore we neglect density-dependent terms in our analysis because they scale with Ma 2 . The coefficients t and t are assumed to be scalars and were obtained from linear fits between time-averaged enth and and between and , respectively. Results from our simulations are shown in figure 1. We normalize t and t by which is an order of magnitude estimate for the turbulent diffusion coefficients. We note that in the parameter regimes studied here, the estimates t0 and t0 are very similar in all of our runs. Our results show that for low Péclet numbers the turbulent heat diffusion increases in proportion to Pe for Pe 1. This is consistent with earlier numerical results for turbulent viscosity (e.g. Käpylä et al. 2020), magnetic diffusivity (e.g. Sur et al. 2008), and passive scalar diffusion (Brandenburg et al. 2009), and with corresponding analytic results in the diffusion dominated (Pe 1) regime. For sufficiently large Pe,˜t tends to a constant value. The turbulent viscosity is also roughly constant in the parameter space covered here. For low fluid Reynolds numbers t is proportional to Re as has been shown in Käpylä et al. (2020). However, we do not cover this parameter regime with the current simulations.

Turbulent viscosity and heat diffusion from decay experiments
Decay experiments were made as an independent way to measure the turbulent viscosity and heat diffusion. Snapshots from the imposed velocity/entropy gradient runs were used as initial conditions and the relaxation terms of the rhs of the Navier-Stokes and entropy equations were deactivated. The large-scale velocity and entropy profiles in such runs decay due to the combined effect of molecular and turbulent diffusion. To measure the decay rate we monitored the amplitude of the = 1 components of and . Exponential decay laws The 1 components of these fields decay exponentially when the forcing is turned off; see panels (c) and (d) of figure 2. Ultimately the amplitude of the 1 mode decreases sufficiently such that it cannot be distinguished from the background turbulence. The time it takes to reach this state varies and depends on the initial amplitudes 0 and 0 . However, at the same time these amplitudes need to be kept as low as possible to avoid non-linear effects becoming important (see, e.g. Käpylä et al. 2020). This is particularly important for the velocity field due to which the range from which turbulent viscosity can be estimated is limited which necessitates running several experiments with different snapshots as initial conditions to reach converged values for t and t .
Due to this, only a limited subset of the parameter range covered by the imposed cases were repeated with decay experiments. We used ten snapshots from each run for the decay experiments. The separation between the snapshots is roughly Δ = 40 rms f such that the realizations can be considered uncorrelated. Results from the decay experiments are shown in figure 1 with crosses ( t ) and pluses ( t ). We find that the results from the decay experiments are consistent with those from the imposed flow and entropy gradient methods. where = 2 /2 is the turbulent kinetic energy, = 2 is the variance of the temperature fluctuations, and are assumed to be universal constants †. Viscous and thermal dissipation rates are defined as = 2 [ S 2 − S 0 2 ] and = (∇ ) 2 = [ (∇ ) 2 − (∇ ) 2 ], respectively, where we have removed contributions from the mean flow and the mean entropy; S 0 denotes the traceless rate-of-strain tensor as defined in (2.4) but with 0 instead of . Pandey et al. (2021) computed Pr t using the − model by fixing the ratio of / , which yields (3.7) For simplicity, we assume / = 1 in this subsection. The results are shown in figure 3.
We find that taking the ratio / to be a constant leads to results where Pr ( − ) t increases monotonically with decreasing Pr when Pe is larger than about 20; see the inset in figure 3 which reveals a dependence of Pr −0.25 . Qualitatively, this result is in agreement with the one in Pandey et al. (2021). However, we would like to note here that a strong assumption was made to reach this conclusion, namely, that the ratio / is fixed, and that and are universal constants, independent of control parameters such as Re and Pe. Henceforth, we relax these assumptions, and also omit primes from the coefficients and . † Primes indicate the constancy of and , but see section 3.4 where this assumption is lifted.

Relaxing the assumption that and are constants
It is reasonable to assume that for sufficiently large Reynolds and Péclet numbers and tend to constant values. Furthermore, there is evidence from numerical simulations that also tends to a non-zero constant value for large Reynolds numbers. Similar evidence for has not been presented. Therefore it is not clear whether the assumption of universality of and is valid. This is particularly important for numerical simulations such as those in the current study where the Reynolds and Péclet numbers are still modest. Since we have independently measured t and t using the imposed flow and entropy method (section 3.1) and from decay experiments (section 3.2), we can estimate and using: where t and t are the ones obtained above in section 3.1 with the imposed field method. The results are shown in figure 4. Our results indicate that and are highly variable and that they depend not only on Re and Pe but also on Pr. Furthermore, for sufficiently large Reynolds and Péclet numbers, and show decreasing trends proportional to roughly −0.25 power of Re and Pe, respectively. This shows that any estimate of t or t with the − model in the parameter regime studied here would require prior knowledge of and for the particular parameters (Re, Pe, Pr) of that system.

Turbulent Prandtl number
Our results for the turbulent Prandtl number Pr t = t t (3.10) are shown in figure 5 with and as discussed in sections 3.1 and 3.2. We find that for Pe 1 the turbulent Prandtl number is roughly inversely proportional to Pe for low molecular Prandtl number. We have not computed the turbulent Prandtl number for cases where both Re and Pe are smaller than unity. For sufficiently high Péclet number the turbulent Prandtl number tends to a constant value which is close to 0.7. This is in accordance with theoretical estimates in that Pr t is somewhat smaller than unity. For example, Rüdiger (1989) derived Pr t = 2/5 using first-order smoothing approximation. The turbulent Prandtl number plays an important role also in the atmospheric boundary layer where several methods yield values of the order of unity (Li 2019, and references therein). That, the turbulent Prandtl number Pr t reaches a constant value at sufficiently large Pe, independent of Pr, is in stark contrast to the results obtained from the − model with a fixed / ; compare figures 3 and 5. Now we make an attempt to understand the reason for this discrepancy. From figure 4 we note the following approximate scaling relations at sufficiently large Re and Pe: ∝ Re −0.25 and ∝ Pe −0.25 , suggesting thus that the ratio / scales with the Prandtl number as Pr +0.25 . With this, if we let / ∝ Pr +0.25 in 3.7, instead of a fixed ratio, and note from the inset of figure 3 that the factor / ∝ Pr −0.25 , we would obtain from 3.7 that Pr ( − ) t becomes independent of Pr, agreeing thus qualitatively with our results as shown in figure 5. Therefore we conclude that the results from the − model with a fixed value for the ratio / are unreliable and that the strong dependence of Pr t on Pr found in Pandey et al. (2021) is due to the restrictive assumption in their model.

Conclusions
Using simulations of weakly compressible isotropically forced turbulence with imposed large-scale gradients of velocity and temperature, and corresponding decay experiments, we find that the turbulent Prandtl number Pr t is roughly 0.7 and independent of the microscopic Prandtl number Pr provided that the Peclét number is higher than about ten. This is in stark contrast from the recent results of Pandey et al. (2021) who found that Pr t ∝ Pr −1 from non-Boussinesq simulations of convection. Although the physical setups are quite different, we were able to qualitatively reproduce their finding under the strong assumption that / is fixed, and note that the method by which the turbulent viscosity and thermal diffusivity were obtained in Pandey et al. (2021) produce unreliable results even in the simpler cases considered here. Relaxing their assumption of universality of and , we find that these depend not only on Re and Pe, respectively, but also on Pr. This allows us to understand the reason for the discrepancy.