Highlight
Free Access
Issue
A&A
Volume 528, April 2011
Article Number A76
Number of page(s) 9
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201014447
Published online 04 March 2011

© ESO, 2011

1. Introduction

Since mixing in stars is a complex interplay of regimes of unstable stratification, stable stratification, differential rotation, gravity waves, double-diffusion, etc., it could be expected that the formalisms employed are general enough to account for such a wide variety of processes. And yet, the literature shows that this is not the case since the two main methodologies being employed are: a) large scale numerical simulations (e.g., Brummell et al. 2002); and b) heuristic models (e.g., Maeder & Meynet 2001; Mathis et al. 2004; Palacios et al. 2003, 2006; Charbonnel & Talon 2005, 2007).

In the first case, the values of some of the stellar parameters, such as the Prandtl number, are widely different from those in stars, and no specific results have yet been produced that can be employed in stellar structure/evolution calculations. The authors of those studies in fact state that their primary goal was to elucidate the intertwined physical processes and not (yet) to provide stellar studies with tools to model the processes of interest. As for the heuristic models, we believe they may have exhausted their fruitfulness since it is extremely difficult, if not almost impossible, to build models capable of embracing the wide variety of processes such as the ones cited at the beginning, ultimately because such processes are not additive.

It is instructive to point out that an analogous situation once existed in geophysics, specifically in modeling atmospheric and oceanic mixing. A few decades ago, however, a shift took place from heuristic, unpredictive approaches in favor of a predictive, prognostic, and flexible tool known as RSM (Reynolds stress models), which are now commonly used. For reasons unclear to the present author, stellar mixing studies have lagged behind geophysical studies, and this paper will thus present a new RSM-based model, as well as the limitations of the heuristic models.

The RSM consist of a set of equations for the turbulent correlations of the velocity, temperature, and concentration fields that are derived directly from the Navier-Stokes, temperature, and concentration equations. The main features can be summarized as follows:

  • a)

    mathematical structure: most of the relevant equations arelinear, algebraic relations and thus pose no particular numericalproblems;

  • b)

    flexibility: one of the key advantages is that adding new processes such as rotation, vorticity, double-diffusion, etc., does not require guesswork since the RSM have a well-defined set of procedural rules,

  • c)

    assessment: the results of the RSM can be assessed before they are used in a stellar context, an important feature that none of the heuristic model enjoys, raising the justified suspicion that such models, tailored to a specific astrophysical setting, have a limited predictive power. The general properties of the RSM have recently been reviewed in (Canuto 2006, 2009).

2. Mean variables

When dealing with mixing in stellar interiors, most of the models concentrate on the temperature field because it enters the convective flux. However, a single field such as the temperature, cannot represent the complexity of stellar mixing processes for which one needs three fields, ui,T,c% subequation 528 0 \begin{eqnarray} u_{i,} T,c \end{eqnarray}(1a)representing velocity, temperature, and a scalar c-field. The velocity field is needed because its mean values correspond to meridional and azimuthal currents, while its gradients give rise to shear (differential rotation) mixing. The scalar c-field can be of two kinds: active and passive. An active c-field is, for example, the mean molecular weight μ\hbox{$\overline \mu $}, and a passive c-field is for example a tracer such as Li7. By active is meant a scalar that affects the density, which in turn affects the velocity field, and thus ultimately alters the turbulent state; in contrast, a passive scalar is transported by the flow without altering it. This can be seen by considering the relations ρ=ρ+ρ,ρ=ρ(αTT+αμμ)% subequation 528 1 \begin{equation} \label{eq1} \rho = \overline \rho +{\rho }',\;\;\;{\rho }'=\overline \rho (-\alpha _T {T}' + \alpha _\mu {\mu }') \end{equation}(1b)where αT,μ =  − (lnρ / ∂T)p,μ,(lnρ / ∂μ)p,T are the expansion coefficients. A passive scalar does not contribute to the rhs of (1b). When we expand the density field in the derivations, we use (1b) in which the c-field is the μ-field. Their different effects on the density mean that the two c-fields have different dynamics, and in what follows we use the c-field notation when the relations apply to both passive and active scalars. When the differentiation becomes necessary, we use the symbol μ\hbox{$\overline \mu$} for an active scalar. The inclusion of the μ\hbox{$\overline \mu$} field is necessary since its gradient is opposite to that of the T-field, a situation that in the ocean gives rise to two well known processes, namely salt-fingers and diffusive convection, the latter being called semi-convection in the astrophysical literature. Each field has an average and a fluctuating component, u=u+u,T=T+T,c=C+c,% subequation 528 2 \begin{equation} \label{eq2} \vec u = \overline{{\vec u}} +u', T = \overline T + {T}', c = \overline C +c', \end{equation}(1c)and since due to turbulence, the correlation among two fluctuating variables is not zero, in addition to having three equations for the mean fields T,C,u\hbox{$\overline T, \overline C, \overline {\vec u}$}, one must derive the equations for second-order correlations, such as heatfluxes:uT,cfluxes:uc,momentumfluxes:uiuj.% subequation 528 3 \begin{equation} \label{eq88} {\mbox{\rm heat fluxes:}}\ \overline{{\vec u'}T'},\ c- \mbox{\rm fluxes:}\ \overline{{ \vec u'}c'},\ \mbox{\rm momentum fluxes:}\ \overline{u'_i u'_j}. \end{equation}(1d)However, it turns out that the equations for the fluxes (1d) involve the correlations of the fields themselves, specifically temp.variance:T2,cvariance:c2,temp.ccorrelation:cT,% subequation 528 4 \begin{eqnarray} \label{eq3} \mbox{\rm temp. variance}\hspace*{-1.5mm} &:& \overline{T'^2}, c{-}{\rm variance} \nonumber\\ &:&\overline{c'^2}, {\rm temp.}-c{-}{\rm correlation:} \overline{c'T'}, \end{eqnarray}(1e)whose dynamic equations must therefore also be included and solved. In what follows, we first derive the dynamic equations for the three mean variables T,C,u\hbox{$\overline T, \overline C, \overline{{\vec u}}$} and then derive those for the variables (1c, d).

2.1. Mean temperature

Since there is no conservation equation for T, we begin with the entropy S-equation from which one then derives an equation for T. Second, we must consider an entropy equation for a mixture of fluids since we want to include different species, for example He and H. With these provisos, the starting S-equation reads as (Landau & Lifshitz 1987) ρTdSdt=xi(qi+ρ􏽥μJi)+ρJi􏽥μxi+σijuixj% subequation 640 0 \begin{equation} \label{eq4} \rho T \frac{{\rm d}S}{{\rm d}t} = - \frac{\partial }{\partial x_i}(q_i + \rho \widetilde{\mu}J_i )+\rho J_i \frac{\partial \widetilde{\mu}}{\partial x_i }+\sigma _{ij} \frac{\partial u_i}{\partial x_j} \end{equation}(2a)where we have defined the following variables: ddt=∂t+uixi,σij=νρ(ui,j+uj,i)23νρδijuk,k,% subequation 640 1 \begin{equation} \label{eq5} \frac{\rm d}{{\rm d}t}=\frac{\partial}{\partial t}+u_i \frac{\partial }{\partial x_i }, \quad \sigma _{ij} = \nu \rho (u_{i,j} +u_{j,i} )-\frac{2}{3}\nu \rho \delta _{ij} u_{k,k} , \end{equation}(2b)where ui is the total velocity field, ν the kinematic viscosity, ρ the density, and σij the viscous stress tensor. Assuming incompressibility requires that the last term in (2b) is zero. Furthermore, the q-flux is defined as follows: qi=FirρhJiρJiκT􏽥μ∂c|p,T,h=􏽥μT􏽥μ∂T|p,c·% subequation 640 2 \begin{equation} \label{eq6} q_i = F_i^r - \rho{h} J_i - \rho J_i \kappa _T \left. {\frac{\partial \widetilde{\mu }}{\partial c}} \right|_{p,T} ,\quad \quad h = \widetilde{\mu }-T\left. {\frac{\partial \widetilde{\mu }}{\partial T}} \right|_{p,c} \cdot \end{equation}(2c)here, Fir\hbox{$F_i^r$} is the radiative flux which, in the absence of the diffusion flux Ji, would be the only term in qi and 􏽥μ\hbox{$\widetilde{\mu}$} is the chemical potential of the mixture. In a two-fluid model, the densities of the two components are ρc,ρ(1 − c), respectively and the concentration equation reads as dρcdt=∂ρJixi,% subequation 640 3 \begin{equation} \label{eq7} \frac{{\rm d} \rho c}{{\rm d}t} = \frac{\partial \rho J_i }{\partial x_i}, \end{equation}(2d)when the c-field corresponds for example to Li7, the corresponding diffusion flux Ji discussed by Chapman & Cowling (1970), is given by Ji=χc(c,i+κTT-1T,i+κpp-1p,i)% subequation 640 4 \begin{equation} \label{eq8} J_i =\chi _c (c_{,i} +\kappa _T T^{-1}T_{,i} +\kappa _p p^{-1}p_{,i}) \end{equation}(2e)with the standard notation a,i = ∂a / xi, and where κT,κp are dimensionless functions. Often, in the literature Ji includes the density and has the opposite sign, and χc is called D. We have also found that different authors employed different notations, for example Eq. (2.1) of Aller & Chapman (1960), Eq. (3) of Schatzman (1969), Eq. (1) of Michaud (1970), and Eq. (1) of Vauclair & Vauclair (1982). As one can observe from Eq. (2e), the diffusion flux is contributed not only by the first term, which is the traditional gradient of the concentration, but also by the temperature and pressure gradients first introduced by Chapman (1917). Next, one employs the thermodynamic relations:

TdSdt=cpdTdtT􏽥μ∂T|p,cdcdt+Tρ-2∂ρ∂T|p,cdpdt·% subequation 814 0 \begin{equation} \label{eq9} T\frac{{\rm d}S}{{\rm d}t} = c_p \frac{{\rm d}T}{{\rm d}t}-T\left. {\frac{\partial \widetilde{\mu }}{\partial T}} \right|_{p,c} \frac{{\rm d}c}{{\rm d}t}+T\rho ^{-2} \left. {\frac{\partial \rho}{\partial T}} \right|_{p,c} \frac{{\rm d}p}{{\rm d}t} \cdot \end{equation}(3a)Making use of Eq. (2d) and taking Tρ-1∂ρ∂T|p,c=1\hbox{$T \rho ^{-1} \left. {\frac{\partial \rho }{\partial T}} \right|_{p,c} =-1$} corresponding to no radiation pressure, Eq. (2a) becomes ρcpdTdt=dpdtFirxi+σijuixj+χc∂cxi∂hxi% subequation 814 1 \begin{equation} \label{eq10} \rho c_p \frac{{\rm d}T}{{\rm d}t}=\frac{{\rm d}p}{{\rm d}t}-\frac{\partial F_i^r }{\partial x_i }+ \sigma _{ij} \frac{\partial u_i }{\partial x_j }+\chi _{\rm c} \frac{\partial c}{\partial x_i}\frac{\partial h}{\partial x_i } \end{equation}(3b)where we have neglected a term that involves the product of two molecular diffusivities and where h=52kTμ1,2,μ1,2=(μ1-1μ2-1)-1% subequation 814 2 \begin{equation} \label{eq11} h =\frac{5}{2}\frac{kT}{\mu _{1,2} },\quad \mu _{1,2} =\left(\mu _1^{-1} -\mu _2^{-1} \right)^{-1} \end{equation}(3c)where the μ are the mean molecular weights of the two substances. We note that if one employs an equation of state for a perfect gas, the dp/dt term inf (3b) can be absorbed into the left hand side which retains the same form but with cv in lieu of cp. The next step is to employ (1b), separate all the fields into a mean and a fluctuating part, carry out the averages, and obtain the dynamic equation for the mean variables. The necessary steps were discussed in a previous work (Canuto 1999), so we limit ourselves to cite the final result that reads as

ρDDt(cpT+K+Ku+G)=p∂txi(Fir+Ficonv+Fike+ρRijuj)+ρχccxihxi% subequation 890 0 \begin{eqnarray} \label{eq12} \overline \rho \frac{\rm D}{{\rm D}t}(c_p T+K+ K_u +G) &=& \frac{\partial \overline p }{\partial t} - \frac{\partial}{\partial x_i }(\overline {F_i^{\rm r}} + F_i^{\rm conv} \nonumber\\ &\quad +& F_i^{\rm ke} + \overline \rho R_{ij} \overline u _j )+\overline \rho \chi _{\rm c} \frac{\partial \overline c }{\partial x_i }\frac{\partial \overline h }{\partial x_i} \end{eqnarray}(4a)where DDt∂t+ujxj·% subequation 890 1 \begin{equation} \label{eq13} \frac{\rm D}{{\rm D}t} \equiv \frac{\partial }{\partial t}+\overline u _j \frac{\partial }{\partial x_j }\cdot \end{equation}(4b)On the lefthand side of (4a), we have defined the following variables: eddykineticenergy:ρK=12ρuiuimeanfieldkineticenergy:Ku=12uiuigravitationalenergy:giρui=ρDGDt·% subequation 890 2 \begin{eqnarray} \label{eq14} \mbox{eddy kinetic energy:}\quad\overline{\rho} K &=& \frac{1}{2} \overline{\rho{u}'_i {u}'_i} \nonumber\\ \mbox{mean field kinetic energy:}\quad K_u &=& \frac{1}{2} \overline{u _i} \,\,\,\overline{u _i}\\ \mbox{gravitational energy:}\quad g_i \overline{\rho} \overline{u_i} &=& \overline{\rho} \frac{{\rm D}G}{{\rm D}t} \cdot \nonumber \end{eqnarray}(4c)On the right side of (4a), we have the variables: radiativeflux:Fir=Kr∂Txiconvectiveflux:Ficonv=cpρTuifluxofeddykineticenergy:Fike=12ρujujuiReynoldsstresses:ρRij=ρujui,% subequation 890 3 \begin{eqnarray} \label{eq15} \mbox{ radiative flux:}\quad F_i^{\rm r} & = & - K_{\rm r} \frac{\partial T}{\partial x_i} \nonumber\\ \mbox{convective flux:}\quad F_i^{\rm conv} &=& c_{\rm p} \overline{\rho T' u'_i} \nonumber\\ \mbox{flux of eddy kinetic energy:}\quad F_i^{\rm ke} &=& \frac{1}{2} \overline{\rho u'_ j u'_j u'_i} \\ \mbox{Reynolds stresses:}\quad\overline{\rho} R_{ij} &=& \overline{ \rho u'_j u'_i}, \nonumber \end{eqnarray}(4d)where Kr = cpρχ and χ, the thermometric conductivity, has the units of cm2 s-1. The diffusion approximation in the first of (4d) is not required by the present formalism and is being suggested as an example. At first sight, Eq. (4a) may be somewhat surprising since one is accustomed to the much simpler equation ρcvDTDt=xi(Fir+Fic+Fike),% subequation 890 4 \begin{equation} \label{eq16} \overline \rho c_v \frac{{\rm D}T}{{\rm D}t} = - \frac{\partial}{\partial x_i}\left(\overline {F_i^{\rm r}} + F_i^{\rm c} +F_i^{\rm ke}\right), \end{equation}(4e)which in the stationary case further reduces to the “flux conservation law”: Fir+Fic+Fike=const.,% subequation 890 5 \begin{equation} \label{eq102} \overline {F_i^{\rm r}} +F_i^{\rm c} + F_i^{ke} ={\rm const.}, \end{equation}(4f)implying the constancy of the radiative + convective + eddy kinetic energy fluxes, the latter being often neglected. Equation (4a) is considerably more general than (4f) for the following reasons:

  • 1)

    the lhs contains the total energy of the system, as indeedexpected, namely, E=cpT+K+Ku+G% subequation 890 6 \begin{equation} \label{eq17} E= c_p T+K+ K_u +G \end{equation}(4g)which is the sum of enthalpy, turbulent kinetic energy, mean flow kinetic energy and gravitational energy, respectively. Which of the four terms in (4g) is the largest depends on the specific physical problem at hand;

  • in the second term of the rhs of (4a), we have four terms Fir+Fic+Fike+ρRijuj% subequation 890 7 \begin{equation} \label{eq18} \overline {F_i^{\rm r}} +F_i^{\rm c} + F_i^{\rm ke} +\overline \rho R_{ij} \overline u _j \end{equation}(4h)which represent the radiative, convective, and eddy kinetic energy fluxes, while the last term is the flux of the Reynolds stresses by the large scale velocity field uj\hbox{$\overline u _j$}. To get a feeling for what the last term means, consider the z-component of (4h) and the closures to be derived in what follows: Fzc=cpραTHpKh(ad),Rzjuj=KmKu∂z.% subequation 890 8 \begin{equation} \label{eq19} F_z^{\rm c} =\frac{c_p \overline \rho}{\alpha _T H_p }K_{\rm h} (\nabla -\nabla _{\rm ad}), \quad R_{zj} \overline u _j =-K_{\rm m} \frac{\partial K_u }{\partial z}. \end{equation}(4i)Here, Kh,m are heat and momentum diffusivities, Ku the mean flow kinetic energy (see the second of 4c) and we have used the relations: uw=Kmu∂z,vw=Kmv∂z·% subequation 890 9 \begin{equation} \label{eq20} \overline {uw} = - K_{\rm m} \frac{\partial \overline u }{\partial z},\quad \overline {vw} = - K_{\rm m} \frac{\partial \overline v }{\partial z}\cdot \end{equation}(4j)Thus, the ratio of the second to last term in (4h) is given by FzcρRzjuj=1σtcpαTHpadzKu% subequation 890 10 \begin{equation} \label{eq21} \frac{F_z^{\rm c}}{\overline \rho R_{zj} \overline u _j }=-\frac{1}{\sigma _t }\frac{c_p}{\alpha _T H_p}\frac{\nabla -\nabla_{\rm ad} }{\partial _z K_u} \end{equation}(4k)where the ratio Kh/Km is the inverse of the turbulent Prandtl number σt, which in the case of stable stratification ∇−∇ad < 0, is an increasing function of the Richardson number as discussed in detail in Paper II. To assign a numerical value to the mean kinetic energy, one needs the meridional currents possibly derived from helio-seismological studies, as discussed in Paper V;

  • an interesting but never accounted for term, is the last one in (4a) which represents the effect of the gradient of the mean c-field on the temperature field.

2.2. Mean c-field

The starting relation is Eqs. (2d, e). Carrying out the mass averaging process (Favre 1969; Canuto 1997) defined as ρc=ρC,ρcui=ρCui+ρuic,% subequation 1099 0 \begin{equation} \label{eq22} \overline{\rho c} = \overline \rho \overline C , \quad \overline{\rho c u_i} = \overline {\rho} \overline C \overline u _i +\overline {\rho {u}'_i {c}'}, \end{equation}(5a)we obtain ρDCDt=xi(ρJic)+xi(ρχcCxi)·% subequation 1099 1 \begin{equation} \label{eq23} \overline \rho \frac{{\rm D} \overline C }{{\rm D}t} = - \frac{\partial}{\partial x_i}(\overline \rho J_i^c) + \frac{\partial}{\partial x_i} \left( \overline \rho \chi _c \frac{\partial \overline C}{\partial x_i} \right)\cdot \end{equation}(5b)The second term in the rhs of (5b) is the ordinary molecular diffusion, while the first term is turbulent c-flux akin to the convective flux, second term in (4d): cflux:ρJic=ρuic.% subequation 1099 2 \begin{equation} c{-}\mbox{\rm flux:}\quad\overline{\rho} J_ i^{\rm c} = \overline{\rho u'_i c'}. \end{equation}(5c)At first sight, it may look as though Eqs. (5b, c) are unrelated to the T-field but that is not the case since, as we show below in Eqs. (12d, e), the turbulent c-flux (5c) entails the turbulent convective flux. Clearly, the term in (5b) with the molecular diffusivity χc only represents passive scalars.

2.3. Mean velocity

We begin with the Navier-Stokes equations for a compressible fluid: ∂ρui∂t+xjρuiuj=∂pxiρgi+σijxjFi,% subequation 1170 0 \begin{equation} \label{eq24} \frac{\partial \rho u_i}{\partial t} + \frac{\partial}{\partial x_j} \rho u_i u_j = - \frac{\partial p}{\partial x_i} - \rho g_i + \frac{\partial \sigma_{ij}}{\partial x_j} \equiv F_i, \end{equation}(6a)where the viscous stress tensor was defined in Eq. (2b). Mass averaging Eq. (6a) via relations analogous to (5a) ρui=ρui,ρui=0,ρuiuj=ρuiuj+ρuiuj% subequation 1170 1 \begin{equation} \label{eq25} \overline {\rho u_i} = \overline \rho\,\, \overline u _i , \overline {\rho {u}'_i } =0, \overline {\rho u_i u_j} = \overline \rho \overline u _i \overline {u_j} + \overline {\rho {u}'_i {u}'_j } \end{equation}(6b)one obtains from (6a): ρui∂t+xjρ(uiuj+Rij)=pxiρgi% subequation 1170 2 \begin{equation} \label{eq26} \frac{\partial \overline \rho \overline u _i }{\partial t}+ \frac{\partial }{\partial x_j }\, \overline \rho (\overline u _i \overline u _j +R_{ij} )=-\frac{\partial \overline p }{\partial x_i } - \overline \rho g_i \end{equation}(6c)where the Reynolds stresses were defined in Eq. (4d). It is important to notice that Eq. (6c) can be rewritten in the more transparent form: ρDuiDt=xi(pδij+ρRij)ρgi.% subequation 1170 3 \begin{equation} \label{eq27} \overline \rho \frac{{\rm D} \overline u _i }{{\rm D}t}=-\frac{\partial }{\partial x_i }(\overline p \delta _{ij} + \overline \rho R_{ij} ) - \overline \rho g_i . \end{equation}(6d)To solve Eq. (6d) one needs the Reynolds stresses, which we study in Sect. 3. In conclusion, the equations for the mean T,C,u\hbox{$\overline T ,\overline C, \overline {\vec u}$} are given by Eqs. (4a), (5b), and (6d). These equations contain the following second-order turbulent correlations: ρTui,ρcuiρuiuj,% subequation 1170 4 \begin{equation} \label{eq28} \overline {\rho {T}'{u}'_i } ,\quad \overline {\rho {c}'{u}'_i } \quad \overline {\rho {u}'_i {u}'_j } , \end{equation}(6e)representing temperature, c-field, and momentum fluxes whose equations we derive next.

3. Reynolds stresses

The derivation of the dynamic equation for the Reynolds stress was presented in an earlier work (Canuto 1997, Eqs. (14)), so we only cite the final result which we simplify by neglecting the terms due to compressibility. We have DRijDt+Dij=(Rikuj,k+Rjkui,k)􏽼0􏽻􏽺0􏽽shear+ρ-2(ρujδik+ρuiδjk)pxk􏽼0􏽻􏽺0􏽽buoyancy×ρ-1πij􏽼0􏽻􏽺0􏽽pressurecorrealtion23εδij􏽼0􏽻􏽺0􏽽dissipation% subequation 1279 0 \begin{eqnarray} \label{eq29} \frac{{\rm D}R_{ij}}{{\rm D}t}+D_{ij} &=& - \underbrace {(R_{\rm ik} \overline u _{j,k} + R_{jk} \overline u _{i,k})}_{\rm shear}\nonumber\\ & \quad + & \underbrace {\overline \rho \,\, ^{-2} \left(\overline {{\rho }'{u}'_j} \delta_{ik} + \overline {{\rho}'{u}'_i} \delta _{jk}\right) \frac{\partial \overline p}{\partial x_k }}_{\rm buoyancy}\nonumber\\ &\quad \times &\underbrace {-\,\overline \rho\,\, ^{-1}\pi _{ij} }_{\rm pressure\,correaltion}-\underbrace {\frac{2}{3}\varepsilon \delta _{ij} }_{\rm dissipation} \end{eqnarray}(7a)where πij=Πij13δijΠkk\hbox{$\pi _{ij} =\Pi _{ij} -\frac{1}{3}\delta _{ij} \Pi _{kk} $} is the traceless pressure correlation tensor, and Dij represents the diffusion of the Reynolds and pressure fluxes, both of which are third-order moments. The trace of (7a) yields the equation for the turbulent kinetic energy, K=12Rii% subequation 1279 1 \begin{equation} \label{eq30} K=\frac{1}{2}\,R_{ii} \end{equation}(7b)which reads as DKDt+D(K)=Rijui,j+gλi(αTJihαμJiμ)ε% subequation 1279 2 \begin{equation} \label{eq31} \frac{{\rm D}K}{{\rm D}t}+{\rm D}(K) = - R_{ij} \overline u _{i,j} + g \lambda _i \left(\alpha _T J_i^h - \alpha _\mu J_i^\mu \right)- \varepsilon \end{equation}(7c)where λi(gρ)-1pxi·% subequation 1279 3 \begin{equation} \label{eq32} \lambda _i \equiv -\,\left(g\overline \rho \right)^{-1}\frac{\partial \overline p }{\partial x_i }\cdot \end{equation}(7d)Since the density field depends on both T and μ-fields, we also have ρ/ρ=αTT+αμμ% subequation 1279 4 \begin{equation} \label{eq33} {\rho }'/\overline \rho = - \alpha _T {T}' + \alpha _\mu {\mu }' \end{equation}(7e)where αT,μ =  −    (lnρ / ∂T)p,μ,(lnρ / ∂μ)p,T are the expansion coefficients. Thus, the combination of heat and μ-fluxes: Jiρ=ρ-1ρui=g-1bui=αTuiTαμuiμ=αTJihαμJiμ,% subequation 1279 5 \begin{eqnarray} \label{eq34} J_i^\rho &=& - \overline \rho \,\, ^{-1}\overline {{\rho }'{u}'_i }\nonumber\\ &=& g^{-1}\overline {{b}'{u}'_i } =\alpha _T \overline {u_i ^\prime {T}'} - \alpha _\mu \overline {u_i ^\prime {\mu }'} \nonumber\\ & =& \alpha _T J_i^h -\alpha _\mu J_i^\mu, \end{eqnarray}(7f)is called the density-buoyancy flux where b=gρ-1ρ\hbox{$b' = - g \overline \rho \,\, ^{-1}{\rho}'$} is the buoyancy field. Thus, the physical interpretation of the rhs of (7c) is as follows. The first term is the source of K due to mean shear; the second term given by (7f) may act as a source or as a sink depending on whether one deals with unstable or stable stratification, as well as the sign of the μ gradient (in the case of salt fingers; the contribution is a source while in the diffusive convective case, semi-convection, the term acts as a sink); finally, the last term, which is always negative, represents the sink of K due to molecular viscosity, and its equation will be discussed in Sect. 6.

Next, we consider the Reynolds stresses which enter the first source term in the rhs of (7c), and we introduce the traceless form bij=Rij13δijRkk,% subequation 1390 0 \begin{equation} \label{eq35} b_{ij} = R_{ij} - \frac{1}{3}\delta _{ij} R_{kk} , \end{equation}(8a)whose dynamic equation can be derived from Eqs. (7a, c) to be DbijDt+Dij(b)=4K3SijΣijZij+Bijπij,% subequation 1390 1 \begin{equation} \label{eq36} \frac{{\rm D}b_{ij}}{{\rm D}t} + {D}_{ij} (b) =-\frac{4K}{3}S_{ij} -\Sigma _{ij} -Z_{ij} +B_{ij} -\pi _{ij}, \end{equation}(8b)where we have defined the following tensors with zero trace: Shear: Sij=12(ui,j+uj,i),Vorticity: Vij=12(ui,juj,i)% subequation 1390 2 \begin{equation} {\rm Shear:}~S_{ij} =\frac{1}{2}(\overline u _{i,j} +\overline u _{j,i} ), {\rm Vorticity:}~V_{ij} =\frac{1}{2}\left(\overline u _{i,j} -\overline u _{j,i} \right) \end{equation}(8c)Zij=bikVjk+bjkVik,Σij=bikSjk+bjk Sik23δijbkmSkm% subequation 1390 3 \begin{equation} \label{eq37} Z_{ij} = b_{ik} V_{jk} +b_{jk} V_{ik} ,\quad \Sigma _{ij} = b_{ik} S_{jk} +b_{jk}~S_{ik} - \frac{2}{3} \delta_{ij} b_{\rm km} S_{\rm km} \end{equation}(8d)Bij=g(λiJjρ+λjJiρ23δijλkJkρ).% subequation 1390 4 \begin{equation} B_{ij}=g\left(\lambda _i J_j^\rho +\lambda _j J_i^\rho -\frac{2}{3}\delta _{ij} \lambda _k J_k^\rho \right). \end{equation}(8e)The final step is the closure of the pressure correlation tensor, a topic that has been widely discussed elsewhere (Canuto 1994, 1997) and need not be repeated here. The final expression of Eq. (8b) is then DbijDt+Dij(b)=5τbij8K15Sij(1α1)Σij(1α2)Zij+β5Bij.% subequation 1390 5 \begin{eqnarray} \label{eq38} \frac{{\rm D} b_{ij}}{{\rm D} t} + {D}_{ij} (b)=& -& \frac{5}{\tau} b_{ij} - \frac{8K}{15} S_{ij} -(1 - \alpha _1) \Sigma_{ij} \nonumber\\ &-& (1-\alpha_2) Z_{ij} + \beta _5 B_{ij} . \end{eqnarray}(8f)Since α1,2 = 0.984,0.568,   β5 = 1 / 2,τpv = 2τ / 5, we neglect the Σij term and round off α2 = 1 / 2, in which case (8f) simplifies to DbijDt+Dij(b)=5τbij8K15Sij12Zij+12Bij,% subequation 1390 6 \begin{equation} \label{eq39} \frac{{\rm D}b_{ij}}{{\rm D}t} + {D}_{ij} (b)=-\frac{5}{\tau}b_{ij} -\frac{8K}{15}S_{ij} -\frac{1}{2}Z_{ij} +\frac{1}{2}B_{ij}, \end{equation}(8g)which, in the stationary and local limit, becomes bij=875Sij110τZij+110τBij.% subequation 1390 7 \begin{equation} \label{eq40} b_{ij} =-\frac{8}{75}K \tau S_{ij} -\frac{1}{10}\tau Z_{ij} +\frac{1}{10}\tau B_{ij}. \end{equation}(8h)Equation (8h) is quite simple: in fact, it is a set of linear algebraic relations whose solution yields the Reynolds stresses bij. We recall that Zij depends on bij itself, and Bij are given in Sect. 4. We also recall that the dynamical time scale is given by τ=2Kε% subequation 1390 8 \begin{equation} \label{eq41} \tau =\frac{2K}{\varepsilon} \end{equation}(8i)whose determination we discuss in Sect. 6.

4. Heat and μ-fluxes

We begin by defining the heat and μ-fluxes:

ρJih=ρuiTρJiμ=ρuiμ.(9a,b)$$ \overline \rho J_i^h =\overline {\rho {u}'_i {T}'} \quad \overline \rho J_i^\mu = \overline {\rho {u}'_i {\mu }'} .\hspace*{4.2cm} (9{\rm a,b}) $$To derive the dynamic equation for (9a), one begins with Eq. (3b) multiplied by ui and adds to it Eq. (6a) multiplied by T. Neglecting the molecular terms, we have ∂ρTui∂t+xiρTuiuj=Aui+FiT,% subequation 1544 0 \begin{equation} \label{eq43} \frac{\partial \rho Tu_i}{\partial t} + \frac{\partial }{\partial x_i} \rho Tu_i u_j = Au_i +F_i T, \end{equation}(9c)where the variable A is defined as cpAdpdtFirxi+σijuixj,Fi∂pxigiρ+σijxj·% subequation 1544 1 \begin{equation} \label{eq44} c_p A \equiv \frac{{\rm d}p}{{\rm d}t}-\frac{\partial F_i^r}{\partial x_i} + \sigma_{ij} \frac{\partial u_i}{\partial x_j}, \quad F_i \equiv - \frac{\partial p}{\partial x_i} - g_i \rho + \frac{\partial \sigma _{ij}}{\partial x_j}\cdot \end{equation}(9d)Next, we average Eqs. (9c, d). The process was carried out in detail elsewhere (Canuto 1997) with the result DJihDt+Df(Jih)=βjRijJjh(Sij+Vij)+gλi(αTT2αμμT)τ-1Jih,% subequation 1544 2 \begin{eqnarray} \label{eq45} \frac{{\rm D}J_i^h}{{\rm D}t} &+& {\rm D}_f \left(J_i^h \right) = \beta _j R_{ij} -J_j^h (S_{ij} + V_{ij})\nonumber\\ & + & g\lambda _i \left(\alpha_T \overline {{T}'^2} - \alpha _\mu \overline {{\mu}'{T}'}\right) - \tau _{p\theta }^{-1} J_i^h , \end{eqnarray}(9e)where βi=∂Txiλigcp·% subequation 1544 3 \begin{equation} \label{eq46} \beta _i = - \frac{\partial T}{\partial x_i} - \lambda _i \frac{g}{c_p}\cdot \end{equation}(9f)It is useful to consider the physical origin of the third term on the rhs of (9e). Consider the fluctuating component of the term – ρg on the rhs of the momentum Eq. (6a). Since to obtain the dynamical equation for the heat flux (9a), one multiplies the T′ equation by ui\hbox{$u'_i$} and the ui\hbox{$u'_i$} equation by T′, the last operation gives rise to the term gρ-1ρT=g(αTT2αμTμ)% subequation 1544 4 \begin{equation} \label{eq47} - g \overline \rho \, ^{-1} \overline {{\rho }'{T}'} = g \left(\alpha _T \overline {{T}'^2} - \alpha _\mu \overline {{T}'{\mu }'}\right) \end{equation}(9g)where we use (7e). Relation (9g) is indeed the third term in the rhs of (9e). The first term in the rhs of (9e) represents the potential energy (temperature variance), but due to the μ-field, there is a second term representing the correlation between the T and μ fields. In the local and stationary case, Eq. (9e) simplifies considerably as it becomes an algebraic equation for the heat flux: (τ-1δij+Sij+Vij)Jjh=βjRij+gλi(αTT2αμTμ).% subequation 1544 5 \begin{equation} \label{eq48} \left(\tau _{p \theta}^{-1} \delta _{ij} + S_{ij} + V_{ij} \right) J_j^h = \beta _j R_{ij} + g \lambda _i \left(\alpha _T \overline {{T}'^2} -\alpha _\mu \overline {{T}'{\mu }'} \right). \end{equation}(9h)The last two terms can be derived to have the following form: T2=τθJihβi,Tμ=τμθ(JiμβiJihμ,i).% subequation 1544 6 \begin{equation} \label{eq49} \overline {{T}'^2} =\tau _\theta J_i^h \beta _i ,\quad \overline {{T}'{\mu}'} =\tau _{\mu \theta } \left(J_i^\mu \beta _i -J_i^h \overline \mu ,_i\right). \end{equation}(9i)Once substituted in (9h), the latter can be rewritten in the simpler form: HeatFluxes:(δij+ηij)Jjh=γijβj,% subequation 1674 0 \begin{equation} \mbox{Heat Fluxes:}\quad \label{eq50} \left(\delta _{ij} + \eta _{ij} \right) J_j^h = \gamma _{ij} \beta _j , \end{equation}(10a)where we have defined the following tensors: γij=π4τ(Rijπ2gαμτλiJjμ)ηij=π4τ[Sij+Vijλi(π5αTβj+π2αμμ,j)].% subequation 1674 1 \begin{eqnarray} \label{eq51} \gamma _{ij} &=& \pi _4 \tau (R_{ij} -\pi _2 g\alpha_\mu^\tau \lambda _i J_j^\mu ) \nonumber\\ \eta _{ij} &=&\pi _4 \tau [S_{ij} +V_{ij} -g\tau \lambda _i (\pi _5 \alpha _T \beta _j +\pi _2 \alpha _\mu \overline \mu _{,j} )]. \end{eqnarray}(10b)In order to homogenize the notation, we have measured the different time scales in units of (8i) and introduced the following dimensionless variables, π1,2,3,4,5τ/τ,τμθ/τ,τμ/τ,τ/τ,τθ/τ,% subequation 1674 2 \begin{equation} \label{eq52} \pi _{1,2,3,4,5} \equiv \tau _{p\mu } /\tau ,\;\;\tau _{\mu \theta } /\tau ,\;\;\tau _\mu /\tau ,\;\;\tau _{p\theta } /\tau ,\;\;\tau _\theta /\tau , \end{equation}(10c)whose determination is an important task that will be dealt with subsequently. Relations (10c) represent decay times due to smallscale dissipation as well as pressure-correlations time scales (hence the labels “” and “”), which in the turbulence closure literature are known as the Rotta’s (1951) terms. Some comments are needed to interpret the physical content of Eqs. (10a,b). The heat flux depends on the following variables: Adiabatic.temp.gradient,shear/vorticityfluxes.% subequation 1674 3 \begin{equation} \mbox{Adiabatic. temp. gradient, shear/vorticity}, \mu - {\rm fluxes}. \end{equation}(10d)The Reynolds stresses Rij are obtained solving Eqs. (8h), which in turn depend on the density flux, that is, heat fluxes (10a, b) and concentration fluxes that we study next. It is of interest to note that the heat flux in (10a) can be rewritten in the more familiar form given by a heat diffusivity times the adiabatic temperature gradient. Indeed, using the Hamilton-Cayley theorem, one can rewrite (10a) as Jih=(Kh)ijβj% subequation 1674 4 \begin{equation} \label{eq53} J_i^h =(K_{\rm h})_{ij} \beta _j \end{equation}(10e)where (Kh)ij is a turbulent heat diffusivity tensor whose explicit form can be computed, but it is doubtful whether in practical computations such a form could be more useful than (10a). μfluxes,ρJiμ=ρuiμ.$$ \mu {-}{\rm fluxes},\quad \overline {\rm {\rho }} {\rm {J}}_{\rm {i}}^{\rm {\mu }} {\rm {\bf = }}\overline {{\rm {\rho {u}'}}_{\rm {\bf i}} {\rm {\bf {\mu }'}}}. $$ We begin with relations (2d, e). Multiplying (2d) by ui, Eq. (6a) by μ and adding the two, we obtain ∂ρμui∂t+xiρμuiuj=Fiμ.% subequation 1787 0 \begin{equation} \label{eq54} \frac{\partial \rho \mu u_i}{\partial t}+\frac{\partial }{\partial x_i }\rho \mu u_i u_j =F_i \mu . \end{equation}(11a)Next, using (5a), we mass average (11a) and obtain an equation analogous to (9c): DJiμDt+Df(Jiμ)=μ,jRijJjμ(Sij+Vij)+gλi(αTTμαμμ2)τ-1Jiμ.% subequation 1787 1 \begin{eqnarray} \label{eq55} \frac{{\rm D}J_i^\mu }{{\rm D}t} + D_f (J_i^\mu) &=& - \overline \mu,_j R_{ij} - J_j^\mu \left(S_{ij} +V_{ij}\right) \nonumber\\ &\quad +& g\lambda _i (\alpha _T \overline {{T}'{\mu }'} - \alpha _\mu \overline {{\mu }'^2} )- \tau _{p\mu }^{-1} J_i^\mu . \end{eqnarray}(11b)Using the second of (9i) and the analog of the first of (9i) μ2=τμJjμμ,j,% subequation 1831 0 \begin{equation} \label{eq56} \overline {{\mu }'^2} = - \tau _\mu J _j^\mu \overline \mu_{,j}, \end{equation}(12a)we obtain, after some algebra, the following equation for the μ-flux: (δij+ξij)Jjμ=dijμ,j% subequation 1831 1 \begin{equation} \label{eq57} (\delta _{ij} + \xi _{ij} )J_j^\mu = - d_{ij} \overline \mu _{,j} \end{equation}(12b)where dij=π1τ(Rij+gαTπ2τλiJjh)ξij=π1τ[Sij+Vijλi(π2αTβj+π3αμμ,j)].% subequation 1831 2 \begin{eqnarray} \label{eq58} d_{ij} &=& \pi _1 \tau \left(R_{ij} + g \alpha _T \pi _2 \tau \lambda _i J_j^h \right) \nonumber \\ \xi _{ij} & = & \pi _1 \tau \left[S_{ij} +V_{ij} - g \tau \lambda _i \left(\pi _2 \alpha _T \beta _j +\pi _3 \alpha _\mu \overline \mu _{,j}\right) \right]. \end{eqnarray}(12c)As discussed before, Eqs. (10a, b) and (12d, e) show that the heat and μ-fluxes depend on each other and both depend on shear and vorticity. In principle, Eqs. (12d, d) can also be written in a form analogous to Eq. (10e).

5. The dissipation-relaxation time scales π′s

Each of the second-order correlations discussed in the previous sections is characterized by a dissipation-correlation time scale that we have normalized with the dynamical time scale (8i) giving rise to (10c), which must now be determined in terms of the large scale features of the model. Since we are considering a physical situation with both mean velocity fields (giving rise to shear mixing), temperature and μ- fields (giving rise to double-diffusion processes), we need two large-scale, dimensionless variables, one to characterize the temperature-velocity fields and the second to characterize the temperature-μ fields, or, more correctly, the relative gradients of these fields. Traditionally, these two variables have been chosen to be the Richardson number Ri and the density ratio Rμ (a term borrowed from oceanography). Ri is defined as Ri=N2Σ2% subequation 1903 0 \begin{equation} \label{eq59} {\rm Ri} = \frac{N^2}{\Sigma ^2} \end{equation}(13a)where Σ2 is the mean shear squared defined in terms of the mean velocities fields, Σ=(2SijSij)1/2,Sij=12(ui,j+uj,i),ui,j=ui/xj% subequation 1903 1 \begin{equation} \label{eq60} \Sigma =(2S_{ij} S_{ij} )^{1/2},\quad S_{ij} =\frac{1}{2}(\overline u _{i,j} +\overline u _{j,i} ),\quad \overline u _{i,j} =\partial \overline u _i /\partial x_j \end{equation}(13b)while the Brunt-Väisäla frequency N is given by N2=gρρ∂z=gαT[T∂z(T∂z)ad]gαμμ∂z,% subequation 1903 2 \begin{equation} \label{eq61} N^2= -\frac{g}{\rho} \frac{\partial \overline \rho }{\partial z}= g\alpha _T \left[\frac{\partial \overline T }{\partial z}-\left(\frac{\partial \overline T }{\partial z}\right)_{\rm ad} \right] - g \alpha _\mu \frac{\partial \overline \mu }{\partial z}, \end{equation}(13c)where αT,c=(lnρ/T)p,c,(lnρ/C)p,T\hbox{$\alpha _{T,c} = - (\partial \ln \rho / \partial \overline T)_{p,c}, (\partial \ln \rho /\partial \overline C )_{p,T} $} are the expansion coefficients. Using the relations αT[T∂z+(T∂z)ad]=Hp-1(ad)αμμ∂z=Hp-1μ,μlnμlnP·% subequation 1903 3 \begin{eqnarray} \label{eq62} \alpha_T \left[ -\frac{\partial \overline T}{\partial z} + \left(\frac{\partial \overline T}{\partial z}\right)_{\rm ad} \right]&=&H_p^{-1} (\nabla -\nabla _{\rm ad} ) \\[2.5mm] \label{eq63} \alpha _\mu \frac{\partial \overline \mu}{\partial z} &=& - H_p^{-1} \nabla_\mu, \quad \nabla _\mu \equiv \frac{\partial \ln \mu }{\partial \ln P}\cdot \end{eqnarray}the density ratio is defined as Rμgαμμ,zgαT(T,zT,zad)=μad·% subequation 1903 4 \begin{equation} \label{eq64} R_\mu \equiv \frac{g \alpha_\mu \overline \mu_{,z}}{g \alpha _T (\overline T _{,z} - \overline T_{,z \rm ad})} = \frac{\nabla_\mu}{\nabla - \nabla _{\rm ad}}\cdot \end{equation}(13f)Equation (13c) then becomes N2=gHp-1(ad)(1Rμ).% subequation 1903 5 \begin{equation} \label{eq65} N^2 = - gH_p^{-1} (\nabla -\nabla _{\rm ad} )(1-R_\mu). \end{equation}(13g)The task is now that of constructing the functions πk(Ri,Rμ).% subequation 1903 6 \begin{equation} \label{eq66} \pi _k ({\rm Ri}, R_\mu). \end{equation}(13h)The determination of such functions was the main goal of Canuto et al. (2008b, 2009) intended to reproduce laboratory data (without shear but with double-diffusion, DD) and oceanic data (which depend on both shear and DD instabilities). However, in a stellar context, one must include a Peclet number dependence because radiative losses by the eddies weaken the strength of turbulence. Adopting the treatment developed elsewhere (Canuto & Dubovikov 1998; Canuto et al. 2008a,b), we now have: π1=π10(1+Ri Rμa+Rμ)-1,π4=π40f(Pe)(1+Ri1+aRμ)-1π2=π20(1+Ri)-1[1+2RiRμ(1+Rμ2)-1],π5=π50g(Pe),π10=π40=(27 Ko3/5)1/2(1+σt-1)-1,π20=1/3,π3,=π30=π50=σt,f(Pe)=bPe(1+bPe)-1,g(Pe)=cPe(1+cPe)-1,a=10,Ko=5/3,4π2b=5(1+σt-1),7π2c=4σt-1% subequation 1903 7 \begin{eqnarray} \label{eq67} \pi _1 &=& \pi _1^0 \left(1+\frac{{\rm Ri}~R_\mu}{a + R_\mu } \right)^{-1},\,\, \pi _4 = \pi_4^0 f({\rm Pe})\left(1 + \frac{\rm Ri}{1 + a R_\mu}\right)^{-1} \nonumber\\ \pi _2 &=& \pi _2^0 (1 + {\rm Ri})^{-1}[1+ 2 {\rm Ri}\,R_\mu (1 + R_\mu ^2)^{-1}], \,\, \pi _5 = \pi _5^0 g({\rm Pe}), \nonumber\\ \pi _1^0 &=& \pi _4^0 = (27~{\rm Ko}^3/5)^{-1/2}(1 + \sigma _t^{-1} )^{-1}, \pi _2^0 = 1/3, \pi _3, \\ &=& \pi _3^0 = \pi _5^0 = \sigma _t, \nonumber\\ f({\rm Pe}) &=& b{\rm Pe}(1+b{\rm Pe})^{-1},\quad \quad g({\rm Pe})=c{\rm Pe}(1 + c {\rm Pe})^{-1},\nonumber\\ a &=& 10,{\rm Ko} = 5/3,4 \pi ^2b = 5 \left(1+\sigma _t^{-1} \right), 7 \pi ^2c = 4 \sigma _t^{-1} \nonumber \end{eqnarray}(13i)where the Peclet number Pe has the following form: Pe=4π2125K2εχ=1.3wχ-1% subequation 1903 8 \begin{equation} \label{eq68} {\rm Pe} = \frac{4 \pi ^2}{125}\frac{K^2}{\varepsilon \chi} = 1.3\, w\, \ell \chi ^{-1} \end{equation}(13j)where χ (cm2 s-1) was defined after Eq. (4d). Here, is a length scale that approaches κz near the boundaries (κ = 0.4 is the von Karman constant) and becomes a constant fraction of the size of the region in the middle of it. This means that  is not a constant and that decreases near the boundaries of a convective region. This behavior, together with that of w, makes Pe very small near the convective-radiative interface, a conclusion that is model independent. We suggest, however, using the first relation in (13j) with the equations for K-ε given below.

6. The K-ε model

The above relations contain two turbulence variables that must be determined, the turbulent kinetic energy K and its rate of dissipation ε. With the inclusion of compressibility effects, the equation for K reads (Canuto 1997, Eqs. (15b) and (35c)): DKDt+D(K)=bijSij􏽼0􏽻􏽺0􏽽Ps+gλiJiρ􏽼0􏽻􏽺0􏽽Pbε.% subequation 2088 0 \begin{equation} \label{eq69} \frac{{\rm D}K}{{\rm D}t} + {\rm D}(K) = \underbrace {-\,b_{ij} S_{ij} }_{P_s} + \underbrace {g\lambda _i J_i^\rho}_{P_b}- \varepsilon . \end{equation}(14a)The diffusion term D(K) is defined as the divergence of the flux of kinetic energy, D(K)=Fikexi,Fike=12q2ui,q2=uiui,% subequation 2088 1 \begin{equation} \label{eq70} {\rm D}(K) = \frac{\partial F_i^{\rm ke}}{\partial x_i },\quad \quad F_i^{\rm ke} = \frac{1}{2}\overline {q^2{u}'_i} , \quad q^2= {u}'_i {u}'_i, \end{equation}(14b)and its determination requires a closure for the third-order moment Fike\hbox{$F_i^{\rm ke}$}. As we discuss in Paper V, such a nonlocal term becomes indispensable in the context of the OV (overshooting) determination at the bottom of the convective zone. We further propose that, rather than using a specific closure, one can turn the problem around and construct a differential equation for Fike\hbox{$F_i^{\rm ke}$} thus avoiding closure problems. The equation for ε, the rate of dissipation of turbulent kinetic energy K, is given by (c1,2 = 2.88, 3.84): 2KεDεDt=12Fikexi+c1Ps+c3Pbc2ε% subequation 2088 2 \begin{equation} \label{eq71} \frac{2K}{\varepsilon }\frac{{\rm D} \varepsilon }{{\rm D}t} = \frac{1}{2}\frac{\partial F_i^{\rm ke}}{\partial x_i} + c_1 P_s +c_3 P_b - c_2 \varepsilon \end{equation}(14c)where we use the form of the flux of ε suggested by Kupka & Mutsham (2007): Fiε=12τ-1Fike·% subequation 2088 3 \begin{equation} \label{eq72} F_i^\varepsilon = -\frac{1}{2}\tau ^{-1}F_i^{\rm ke}\cdot \end{equation}(14d)

7. Features of the model: no Ri(cr)

Stably stratified flows are characterized by two competing factors: a shear instability that acts as a source of mixing and a stable temperature gradient that acts as a sink. The two combine into a single parameter, the Richardson number Ri defined in Eq. (13a). The key physical question is whether there is a critical value of Ri beyond which stable stratification overpowers the mixing action of shear so as to lead to zero mixing. The standard RSM models used thus far exhibit a finite critical Richardson number Ri(cr). For example, the Mellor & Yamada model (1982) yields Ri(cr) = 0.19, while the recent model by Cheng et al. (2002) increases the value to Ri(cr) = O(1), thus overcoming the difficulty first raised by Martin (1985) that the MY model yielded too shallow ocean mixed layers. A stability analysis including nonlinearities (Abarbanel et al. 1984) also arrived at Ri(cr) = O(1).

However, even an Ri(cr)  ≈  1 is not satisfactory. In fact, several recent data have shown that mixing persists even after Ri = 1. These data include meteorological observations (Kondo et al. 1978; Bertin et al. 1997; Mahrt & Vickers 2005; Uttal et al. 2002; Poulos et al. 2002; Banta et al. 2002), lab experiments (Strang & Fernando 2001; Rehmann & Koseff 2004; Ohya 2001), LES (Zilitinkevich et al. 2007; Zilitinkevich & Esau 2007), DNS (Stretch et al. 2001), oceanic measurements (Mack & Schoeberlein 2004), and theoretical modeling (Sukoriansky et al. 2005).

Zilitinkevich et al. (2007) constructed a model called “energy and flux-budget turbulence closure” that has no Ri(cr). They then show how the model can explain a wide variety of data. Their model differs significantly from the “main stream” closure models used in traditional second-order closure models. In particular, their parameterization of the pressure-correlations differs from the “standard” ones developed by numerous authors and supported by theoretical and experimental data (Rotta 1951; Launder et al. 1975; Lumley 1978; Zeman & Lumley 1979; André et al. 1978; Mellor & Yamada 1982; Cheng et al. 2002).

Given this situation, Canuto et al. (2008a,b) tried to answer the following question: can the traditional RSM models encompass an arbitrarily large Ri(cr) and what changes are needed to do so without spoiling their documented performance for small and medium Richardson numbers? The answer turned out to be relatively simple since a minimum alteration was required to the previous mixing models to encompass a no-Ri(cr) feature. Specifically, the introduction of the 1 + Ri term in π4 given in Eq. (13i) was all that was needed (such correction does not exist in the unstably stratified case however). Several arguments justifying such an extension were presented in Canuto et al. (2008a), where the performance of the model vs. data is fully discussed.

8. Features of the model: arbitrary Peclet number

In the numerical simulations of Brummel et al. (2002) and Brun & Toomre (2002), Pe is defined as in (13j), but the length scale  is taken to be constant (the depth of the domain), which as we have discussed, is not true near the boundaries. This makes it difficult to compare the results of the numerical simulations in which Pe is taken tobe larger then unity, the smallest value being about 10. As an example, it is difficult to physically interpret Fig. 14 in Brummel et al. (2002) showing that the OV extent decreases as Pe increases. With a constant length scale, the Pe used in such calculations depends only on the rms velocity, and a large Pe corresponds to a large rsm that should entail large rather small OV’s. In the present treatment, Pe is a dynamical variable computed along with all the other variables and not a quantity whose value is assumed.

9. Summary

In many ways, the most surprising aspect of the RSM results is the relative simplicity of the equations determining the Reynolds stresses, heat, and concentration fluxes, since they are linear algebraic equations. This is especially relevant if one considers the amount of information they contain: stable stratification, unstable stratification, rigid rotation, shear, and radiative losses (Peclet number). The two differential equations for K and ε characterize any turbulence model, not this one in particular. We derived the equations that describe mixing in the presence of temperature,μ,velocity% subequation 2495 0 \begin{equation} {\rm temperature}, \mu, {\rm velocity} \end{equation}(15a)Summarizing our results, we have derived the following equations: meantemperature,Eq. (4a)$$ \mbox{mean temperature, Eq.~(4a)} $$meanμ,Eq. (5b)$$ \mbox{mean} \, \mu, {\rm Eq.~(5b)} $$meanvelocityfield,Eq. ()$$ \mbox{mean \, velocity \, field, Eq.~(\ref{eq27})} $$Reynoldsstresses,Eq. (8h)% subequation 2495 1 \begin{equation} \mbox{Reynolds stresses, Eq.~(8h)} \end{equation}(15b)heatfluxes,Eqs. (10a,b)$$ {\rm heat\,\,\, fluxes, Eqs.~(10a, b)} $$μfluxes,Eqs. (12d,e)$$ \mu- {\rm fluxes, \, Eqs.~(12d, e)} $$timescales:Eqs. (13i)$$ \mbox{\rm time scales:}\ \rm Eqs.~(13i) $$Pecletnumber,Eq. (13j)$$ {\rm Peclet \, number, Eq.~(13j)} $$Kε,Eqs. (14ad).$$ K-\varepsilon, {\rm Eqs.~(14a{-}d)}. $$The model is valid for arbitrary temperature, concentration, and mean velocity gradients, as well as arbitrary Ri and Peclet numbers.

References

  1. Abarbanel, H. D. I., Holm, D. D., Marsden, J. E., & Ratiu, T. 1984, Phys. Rev. Lett., 52, 2352 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aller, L. H., & Chapman, S. 1960, ApJ, 132, 1 [Google Scholar]
  3. André, J. C., DeMoor, G., Lacarrére, P., Therry, G., & du Vachat, R. 1978, J. Atmos. Sci., 35, 1861 [Google Scholar]
  4. Banta, R. M., Newsom, R. K., Lundquist, J. K., et al. 2002, Boundary-Layer Meteorol., 105, 221 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bertin, F., Barat, J., & Wilson, R. 1997, Radio Science, 32, 791 [NASA ADS] [CrossRef] [Google Scholar]
  6. Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825 [Google Scholar]
  7. Brun, A. S., & Toomre, J. 2002, ApJ, 570, 865 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Canuto, V. M. 1997, ApJ, 482, 827 [Google Scholar]
  9. Canuto, V. M. 1999, ApJ, 524, 311 (cited as C99) [NASA ADS] [CrossRef] [Google Scholar]
  10. Canuto, V. M. 2006, Theoretical modeling of Convection. I. Key physical processes, in Convection in Astrophysics, ed. F. Kupka, I. W. Roxburgh, & K. L. V. Chan, 3, Part II, Reynolds Stress Model, IAU Symp., 239, 19 [Google Scholar]
  11. Canuto, V. M. 2009, Turbulence in astrophysical and geophysical flows, in Interdisciplinary Aspects of Turbulence, Lect. Notes Phys., ed. W. Hillebrandt, & F. Kupka, 107 (Berlin: Springer Verlag) [Google Scholar]
  12. Canuto, V. M., & Dubovikov, M. S. 1998, ApJ, 493, 834 [Google Scholar]
  13. Canuto, V. M., Cheng, Y., Howard, A. M., & Esau, I. N. 2008a, J. Atmos. Sci., 65, 2437 [NASA ADS] [CrossRef] [Google Scholar]
  14. Canuto, V. M., Cheng, Y., & Howard, A. M. 2008b, GRL, 35, L02613 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chapman, S. 1917, MNRAS, 77, 539 [NASA ADS] [Google Scholar]
  16. Chapman, S., & Cowling, T. G. 1970, The mathematical theory of non-uniform gases (Cambridge, Univ. Press) [Google Scholar]
  17. Charbonnel, C., & Talon, S. 2005, Science, 309, 2189 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  18. Charbonnel, C., & Talon, S. 2007, Science, 318, 922 [CrossRef] [PubMed] [Google Scholar]
  19. Charbonnel, C., & Zahn, J. P. 2007, A&A, 467, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Cheng, Y., Canuto, V. M., & Howard, A. M. 2002, J. Atmos. Sci., 59, 1550 [NASA ADS] [CrossRef] [Google Scholar]
  21. Favre, A. 1969, in Problems of Hydrodynamics and continuous mechanics, Philadelphia, SIAM, 231 [Google Scholar]
  22. Kondo, J., Kanechika, O., & Yasuda, N. 1978, J. Atmos. Sci., 35, 1012 [NASA ADS] [CrossRef] [Google Scholar]
  23. Kupka, F., & Muthsam, H. J. 2007, in Convection in Astrophysics, ed. F. Kupka, I. W. Roxburgh, & K. L. Chan, IAU Symp., 239, 86 (Cambridge Univ. Press) [Google Scholar]
  24. Landau, L. D., & Lifshitz, E. M. 1987, Fluid Mechanics (New York: Pergamon Press) [Google Scholar]
  25. Launder, B. E., Reece, G., & Rodi, W. 1975, J. Fluid Mech., 68, 537 [NASA ADS] [CrossRef] [Google Scholar]
  26. Lumley, J. L. 1978, Adv. Appl. Mech., 18, 123 [Google Scholar]
  27. Mack, S. A., & Schoeberlein, H. C. 2004, J. Phys. Oceanogr., 34, 736 [NASA ADS] [CrossRef] [Google Scholar]
  28. Maeder, A., & Meynet, G. 2001, A&A, 373, 555 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Mahrt, L., & Vickers, D. 2005, Boundary-Layer Meteorol., 116, 313 [NASA ADS] [CrossRef] [Google Scholar]
  30. Martin, P. J. 1985, J. Geophys. Res., 90, 903 [NASA ADS] [CrossRef] [Google Scholar]
  31. Mathis, S., Palacios, A., & Zahn, J. P. 2004, A&A, 425, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Mellor, G. L., & Yamada, T. 1982, Rev. Geophys. Space Phys., 20, 851 [NASA ADS] [CrossRef] [Google Scholar]
  33. Michaud, G. 1970, ApJ, 160, 641 [NASA ADS] [CrossRef] [Google Scholar]
  34. Ohya, Y. 2001, Boundary-Layer Meteorol., 98, 57 [NASA ADS] [CrossRef] [Google Scholar]
  35. Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Palacios, A., Charbonnel, C., Talon, S., & Siess, L. 2006, A&A, 453, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Poulos, G. S., Blumen, W., Fritts, D. C., et al. 2002, Bull. Amer. Meteor. Soc., 83, 555 [CrossRef] [Google Scholar]
  38. Rehmann, C. R., & Koseff, J. R. 2004, Dynamics of Atmospheres and Oceans, 37, 271 [NASA ADS] [CrossRef] [Google Scholar]
  39. Rotta, J. C. 1951, Z. Physik, 129, 547 [Google Scholar]
  40. Schatzman, E. 1969, A&A, 3, 331 [NASA ADS] [Google Scholar]
  41. Strang, E. J., & Fernando, H. J. S. 2001, J. Phys. Ocean., 31, 2026 [Google Scholar]
  42. Stretch, D. D., Rottman, J. W., Nomura, K. K., & Venayagamoorthy, S. K. 2001, 14th Australasian Fluid Mechanics Conference, Adelaide, Australia, 10–14 December [Google Scholar]
  43. Sukoriansky, S. B., & Galperin, I. 2005, Phys. Fluids, 17, 1 [CrossRef] [Google Scholar]
  44. Uttal, T., Curry, J. A., McPhee, M. G., et al. 2002, Bull. American Meteorol. Soc., 83, 255 [NASA ADS] [CrossRef] [Google Scholar]
  45. Vauclair, S., & Vauclair, G. 1982, ARA&A, 30, 3 [Google Scholar]
  46. Zeman, O., & Lumley, J. L. 1979, Turbulent Shear Flows (New York: Springer), 1, 295 [Google Scholar]
  47. Zilitinkevich, S., & Esau, I. 2007, Boundary-Layer Meteorol., 125, 193 [NASA ADS] [CrossRef] [Google Scholar]
  48. Zilitinkevich, S. S., Elperin, T., Kleeorin, N., & Rogachevskii, I. 2007, Boundary Layer Meteorol., 125, 167 [Google Scholar]

Appendix A: Summary of the model equations

Reynolds stresses: bij=Rij13δijRkk,ρRij=ρuiujbij=875Sij110τZij+110τBij\appendix \setcounter{section}{1} \begin{eqnarray} \label{eq73} b_{ij} &=& R_{ij} - \frac{1}{3}\delta _{ij} R_{kk} ,\quad \quad \overline \rho R_{ij} =\overline {\rho {u}'_i {u}'_j } \\ \label{eq74} b_{ij} &=& -\frac{8}{75}K\tau S_{ij} -\frac{1}{10}\tau Z_{ij} +\frac{1}{10}\tau B_{ij} \end{eqnarray}Sij=12(ui,j+uj,i),Vij=12(ui,juj,i),Zij=bikVjk+bjkVik,BijBijhBijμ\appendix \setcounter{section}{1} \begin{eqnarray} \label{eq75} S_{ij} &=&\frac{1}{2}(\overline u _{i,j} +\overline u _{j,i} ),\quad V_{ij} =\frac{1}{2}(\overline u _{i,j} -\overline u _{j,i} ), \nonumber\\ Z_{ij} &=&b_{ik} V_{jk} +b_{jk} V_{ik} , \quad B_{ij} \equiv B_{ij}^h - B_{ij}^\mu \end{eqnarray}(A.3)Bijh=gαT(λiJjh+λjJih23δijλkJkh),Bijμ=gαμ(λiJjμ+λjJiμ23δijλkJkμ)\appendix \setcounter{section}{1} \begin{eqnarray} \label{eq76} B_{ij}^h &=& g\alpha _T (\lambda _i J_j^h +\lambda _j J_i^h -\frac{2}{3} \delta _{ij} \lambda _k J_k^h ),\nonumber\\ B_{ij}^\mu &=& g\alpha _\mu \left(\lambda _i J_j^\mu +\lambda _j J_i^\mu -\frac{2}{3}\delta _{ij} \lambda _k J_k^\mu \right) \end{eqnarray}(A.4)λi(gρ)-1pxi,τ=2Kε·\appendix \setcounter{section}{1} \begin{equation} \label{eq77} \lambda_i \equiv - (g \overline \rho)^{-1} \frac{\partial \overline p}{\partial x_i}, \quad \tau = \frac{2K}{\varepsilon} \cdot \end{equation}(A.5)Heat fluxes,ρJih=ρuiT\appendix \setcounter{section}{1} \hbox{$\overline{\rho} J_{\rm i}^{h} = \overline{\rho {u}'}_{\rm i} {T'}$}: (δij+ηij)Jjh=γijβj,βi=∂Txiλigcpγij=π4τ(Rijπ2gαμτλiJjμ)ηij=π4τ[Sij+Vijλi(π5αTβj+π2αμμ,j)].\appendix \setcounter{section}{1} \begin{eqnarray} \label{eq78} (\delta _{ij} + \eta _{ij})J_j^h &=& \gamma _{ij} \beta _j, \quad \beta _i =-\frac{\partial T}{\partial x_i }-\lambda _i \frac{g}{c_p }\\ \label{eq79} \gamma _{ij} &=& \pi _4 \tau (R_{ij} -\pi _2 g\alpha _\mu \tau \lambda _i J_j^\mu ) \nonumber\\ \eta _{ij} &=& \pi _4 \tau \left[S_{ij} +V_{ij} -g\tau \lambda _i (\pi _5 \alpha _T \beta _j +\pi _2 \alpha _\mu \overline \mu _{,j} )\right]. \end{eqnarray}μ − fluxes, ρJiμ=ρuiμ\appendix \setcounter{section}{1} \hbox{$\overline {\rho} {J}_{i}^{\mu} = \overline {\rho {u}'}_{\rm i} {\mu }'$}: (δij+ξij)Jjμ=dijμ,j\appendix \setcounter{section}{1} \begin{equation} \label{eq80} \left(\delta _{ij} +\xi _{ij} \right) J_j^\mu = - d_{ij} \overline \mu _{,j} \end{equation}(A.8)where: dij=π1τ(Rij+gαTπ2τλiJjh)ξij=π1τ[Sij+Vijλi(π2αTβj+π3αμμ,j)].\appendix \setcounter{section}{1} \begin{eqnarray} \label{eq81} d_{ij} &=& \pi _1 \tau (R_{ij} +g\alpha _T \pi _2 \tau \lambda _i J_j^h ) \nonumber\\ \xi _{ij} &=& \pi _1 \tau [S_{ij} +V_{ij} -g\tau \lambda _i (\pi _2 \alpha _T \beta _j +\pi _3 \alpha _\mu \overline \mu _{,j} )]. \end{eqnarray}(A.9)

Appendix B: Mixing length theory

In the absence of mean velocities and concentration fields, we write the heat flux in the standard form: wT=βχΦ,Φ=Khχ\appendix \setcounter{section}{2} \begin{eqnarray} \label{eq82} \overline {{w}'{T}'} = \beta \chi \Phi , \quad \Phi =\frac{K_h }{\chi} \end{eqnarray}(B.1)where the dimensionless Φ represents the ratio of the turbulent diffusivity to the thermometric one. From Eq. (10a) in the 1D case with j = 3, we derive the relation: wT=Khβ,Kh=τw2π4(1+π4π5Nh2τ2)-1,Nh2=gαTβ.  \appendix \setcounter{section}{2} \begin{equation} \label{eq83} \overline {{w}'{T}'} =K_h \beta , \ K_h =\tau \overline {w^2} \pi _4 \left(1+\pi _4 \pi _5 N_h^2 \tau ^2\right)^{-1},\ N_h^2 = - g\alpha _T \beta.~~ \end{equation}(B.2)Using the relations: w2K~const.,ε=gαTwT,ε=K3/2Λ\appendix \setcounter{section}{2} \begin{equation} \label{eq84} \frac{\overline {w^2} }{K} \sim {\rm const.},\quad \varepsilon =g\alpha_T \overline {{w}'^{T'}} ,\quad \varepsilon =\frac{K^{3/2}}{\Lambda } \end{equation}(B.3)we obtain after some simple algebra (and taking all the coefficients to be of order unity): wT=χβS1/2,Φ~S1/2,S=gαTβΛ4χ2\appendix \setcounter{section}{2} \begin{equation} \label{eq85} \overline {{w}'{T'}} =\chi \beta S^{1/2}, \quad \Phi \sim S^{1/2},\, S=\frac{g\alpha_T \beta \Lambda ^4}{\chi ^2} \end{equation}(B.4)which is the well known expression for the heat flux provided by the MLT for the case of efficient convection. In the same limit, we have that the Peclet number is given by: Pe~(SΦ)1/3~S1/2\appendix \setcounter{section}{2} \begin{equation} \label{eq77} Pe \sim (S\Phi )^{1/3} \sim S^{1/2} \end{equation}(B.5)

Appendix C: The Ri  →  ∞ case, analytic solution

Though this version of the model is restrictive, it is of interest to have a complete analytical solution of the problem which we present in what follows. From relations (10a, b) and (12d, e), we obtain that the heat and μ-fluxes are given by the following expressions: Jh=ρwTρ=Khβ,Kh=τw2AhJμ=ρwμρ=Kμμ∂z,Kμ=τw2Aμ\appendix \setcounter{section}{3} \begin{eqnarray} \label{eq86} J_h &=& \frac{\overline {\rho {w}'{T}'}}{\overline \rho }=K_h \beta , \quad K_h =\tau \overline {w^2} A_h \nonumber\\[-1.5mm] \\ J_\mu &=& \frac{\overline {\rho {w}'{\mu }'}}{\overline \rho} = - K_\mu \frac{\partial \overline \mu}{\partial z},\quad K_\mu =\tau \overline {w^2} A_\mu\nonumber \end{eqnarray}(C.1)where the dimensionless functions Ah,μ are given by the algebraic relations: DAh=π4(1+ηx+π1π2xRμ),DAμ=π1(1+μxπ2π4x),D=(1+ηx)(1+μx)+π1π22π4x2 Rμ,η=π1(π2π3 Rμ),μ=π4(π5π2 Rμ).\appendix \setcounter{section}{3} \begin{eqnarray} \label{eq87} DA_{\rm h} &=& \pi _4 (1 + \eta x + \pi _1 \pi _2 x R_\mu),\nonumber\\ DA_\mu & = & \pi _1 (1 + \mu x - \pi _2 \pi _4 x),\nonumber\\ D &=& (1 + \eta x) (1 + \mu x) + \pi _1 \pi _2^2 \pi _4 x^2~R_\mu, \nonumber\\ \eta &=& \pi _1 (\pi _2 - \pi _3~R_\mu), \nonumber\\ \mu &=&\pi _4 (\pi _5 -\pi _2~R_\mu). \end{eqnarray}(C.2)The variables πs are given by relations (13i) and x is defined as follows: x=τ2Nh2,Nh2=gαTβ=gHp-1(ad).\appendix \setcounter{section}{3} \begin{equation} \label{eq88} x = \tau^2 N_{\rm h}^2 ,\quad N_h^2 = - g \alpha _T \beta = - g H_{\rm p}^{-1} (\nabla -\nabla _{\rm ad}). \end{equation}(C.3)The next step is to solve the Reynolds stress Eqs. (8h). The result is the algebraic relation: w22K=[307+x(AhAμRμ)]-1·\appendix \setcounter{section}{3} \begin{equation} \label{eq89} \overline {\frac{w^2}{2K}} =\left[\frac{30}{7}+ x(A_h -A_\mu R_\mu )\right]^{-1}\cdot \end{equation}(C.4)The two key variables K-ε still remain to be determined and in principle one has to solve the two differential Eqs. (14a, c). The time scale τ then follows from the definition relation (8i) and all the variables are then known as a function of Ri and Rμ.

If a local model can be used, a simplification occurs since Eq. (14a) becomes, after using (7f) and the above results, the following simple expression: ε=τw2(AμRμAh)Nh2.\appendix \setcounter{section}{3} \begin{equation} \label{eq90} \varepsilon = \tau \overline {w^2} (A_\mu R_\mu -A_{\rm h}) N_{\rm h}^2 . \end{equation}(C.5)Combining (1e) with (C.5), one obtains an algebraic relation only in terms of the variable x which reads: x(AμRμAh)=157\appendix \setcounter{section}{3} \begin{equation} \label{eq91} x(A_\mu R_\mu -A_{\rm h}) =\frac{15}{7} \end{equation}(C.6)which, after using relations (1b), can be solved to provide the function: x=x(Ri,Rμ).\appendix \setcounter{section}{3} \begin{equation} x=x(Ri, R_\mu). \end{equation}(C.7)The explicit form of relation (C.6) upon use of relations (1b), reads as follows: x2A(x)+xB(x)157=0\appendix \setcounter{section}{3} \begin{equation} \label{eq92} x^2A(x) + xB(x)-\frac{15}{7}=0 \end{equation}(C.8)where: A(x)=π1(μπ2π4)Rμπ4(η+π1π2Rμ)157(μη+π1π22π4Rμ)B(x)=π1Rμπ4157(η+μ).\appendix \setcounter{section}{3} \begin{eqnarray} A(x) &=& \pi _1 (\mu -\pi _2 \pi _4 )R_\mu - \pi _4 (\eta +\pi _1 \pi _2 R_\mu)\nonumber\\ &\quad -&\frac{15}{7}(\mu \eta +\pi _1 \pi _2^2 \pi _4 R_\mu )\nonumber\\ \label{eq93} B(x)&=& \pi _1 R_\mu - \pi _4 -\frac{15}{7}(\eta +\mu ). \end{eqnarray}(C.9)With the knowledge of x or equivalently τ = 2K / ε, we know the ratio of the variables K-ε but not them singularly. To go further, we adopt a relation of the Kolmogorv type: ε=K3/2Λ\appendix \setcounter{section}{3} \begin{eqnarray} \label{eq94} \varepsilon =\frac{K^{3/2}}{\Lambda } \end{eqnarray}(C.10)where Λ is a mixing length on which we shall comment below. If so, we obtain that the kinetic energy K can now be expressed in terms of the variable x, see relation (C.3), as follows: K=4Λ2(Nh2x)·\appendix \setcounter{section}{3} \begin{equation} \label{eq95} K=4\Lambda ^2 \left(\frac{N_h^2 }{x}\right)\cdot \end{equation}(C.11)If we substitute relation (C.6) into (C.4), we obtain the simple relation: τw2=2815K2ε\appendix \setcounter{section}{3} \begin{equation} \label{eq96} \tau \overline {w^2} = \frac{28}{15}\frac{K^2}{\varepsilon } \end{equation}(C.12)and thus the expressions for the heat and μ diffusivities defined in (C.1) take the forms: Kh=5615Λ2Ah(Nh2x)1/2Kμ=5615Λ2Aμ(Nh2x)1/2·\appendix \setcounter{section}{3} \begin{eqnarray} \label{eq97} K_h &=& \frac{56}{15} \Lambda ^2A_h \left(\frac{N_h^2 }{x}\right)^{1/2} \nonumber\\ K_\mu & = & \frac{56}{15}\Lambda ^2A_\mu \left(\frac{N_h^2 }{x}\right)^{1/2}\cdot \end{eqnarray}(C.13)Relations (C.13) represent the analytical solution of the turbulent part of the problem. See C99, Figs. 1–16.

The last part of the problem that remains to be solved is the determination of the temperature gradient in the presence of turbulent fluxes which is obtained by solving the flux equation which in simplest form is given by: Fr+Fh=const.\appendix \setcounter{section}{3} % subequation 2996 0 \begin{equation} F^{r} + F ^{h} ={\rm const.} \end{equation}(C.14a)where we employ the relations: Fr=Hp-1KrT,Kr=cpρχ,Fh=cpρβKh.\appendix \setcounter{section}{3} % subequation 2996 1 \begin{equation} \label{eq98} F^{\rm r} = H_p^{-1} K_r T\nabla ,\quad K_{\rm r} = c_p \rho \chi , \quad F^{\rm h} = c_{\rm p} \rho \beta K_{\rm h}. \end{equation}(C.14b)Semi-convectionAs discussed in detail in Paper II, in this case we have ∇ − ∇ad > 0,∇μ > 0,Rμ > 0, Nh2<0\appendix \setcounter{section}{3} \hbox{$N_{\rm h}^2 < 0$}. Because of the last relation, the solution of Eq. (C.8) must be taken to be negative x < 0 so that relations (C.13) yield a positive diffusivity. Next, if we introduce the dimensionless variable: U=(adrad)1/2\appendix \setcounter{section}{3} % subequation 2996 2 \begin{equation} \label{eq99} U = \left(\frac{\nabla - \nabla _{\rm ad}}{\nabla _r -\nabla _{\rm ad}}\right)^{1/2} \end{equation}(C.14c)the flux conservation law (C.14a) acquires the form: U2(1+Khχ)=1\appendix \setcounter{section}{3} % subequation 2996 3 \begin{equation} \label{eq100} U^2\left(1+\frac{K_{\rm h}}{\chi}\right)=1 \end{equation}(C.14d)where the ratio Khχ\appendix \setcounter{section}{3} \hbox{$\frac{K_{\rm h}}{\chi}$} measures the importance of turbulent diffusivity over the radiative one (we recall that χ has dimensions of cm2 s-1). Using the first of (C.13), we further have that the ratio Kh / χ can be rewritten as: Khχ=1753π2Ah(x)(x)1/2ΓU,Γ=8π2125[gΛ4χ2Hp(rad)]1/2\appendix \setcounter{section}{3} % subequation 2996 4 \begin{eqnarray} \label{eq101} \frac{K_{\rm h}}{\chi} &=& \frac{175}{3\pi ^2}A_h (x)(-x)^{-1/2}\Gamma U,\nonumber\\ \Gamma &=& \frac{8 \pi ^2}{125}\left[\frac{g\Lambda ^4}{\chi ^2H_{\rm p}}(\nabla _r -\nabla _{\rm ad} )\right]^{1/2} \end{eqnarray}(C.14e)where Γ can be interpreted as an efficiency factor similar to the one in the mixing length scheme (see Appendix B), which can be computed off line. Substituting (C.14e) into (C.14d) yields the following equation for the variable of interest U: ψU3+U2=1,ψ1753π2Ah(x)(x)1/2Γ.\appendix \setcounter{section}{3} % subequation 2996 5 \begin{eqnarray} \label{eq102} \psi U^3+U^2 = 1, \quad \psi \equiv \frac{175}{3\pi ^2}A_{\rm h} (x)(-x)^{-1/2}\Gamma. \end{eqnarray}(C.14f)Finally, we note that the Peclet number is defined in Eq. (13j) and enters in the definitions of the πs given by Eqs. (13i). It is also a function of the variable U since (13j) becomes: Pe=Γ(x)1/2U.\appendix \setcounter{section}{3} % subequation 2996 6 \begin{equation} \label{eq103} {\rm Pe} = \Gamma (- x)^{1/2}U. \end{equation}(C.14g)Salt FingersIn this case we have (see Paper II) ∇ − ∇ad < 0,   ∇μ < 0,   Rμ > 0,Nh2>0.\appendix \setcounter{section}{3} \hbox{$N_{\rm h}^2 > 0.$} Thus, x > 0 and relations (C.11), (C.14c) and  (C.14g) become: K=4gΛ2Hp-1x-1(ad),U=(adadr)1/2,Pe=ΓUx1/2.\appendix \setcounter{section}{3} % subequation 2996 7 \begin{eqnarray} \label{eq104} K &=& 4g\Lambda ^2H_p^{-1} x^{-1}(\nabla _{\rm ad} -\nabla ),\nonumber\\ U &=& \left(\frac{\nabla _{\rm ad} -\nabla }{\nabla _{\rm ad} -\nabla _r }\right)^{1/2},\quad {\rm Pe} = \Gamma Ux^{-1/2}. \end{eqnarray}(C.14h)We further have that: Γ=8π2125[gΛ4χ2Hp(adr)]1/2.\appendix \setcounter{section}{3} % subequation 2996 8 \begin{equation} \label{eq105} \Gamma =\frac{8\pi ^2}{125}\left[\frac{g\Lambda ^4}{\chi ^2H_p }(\nabla _{\rm ad} -\nabla _r )\right]^{1/2}. \end{equation}(C.14i)The first of Eq. (C.14f) does not change but the variable ψ is now defined as follows: ψ1753π2Ah(x)x1/2Γ.\appendix \setcounter{section}{3} % subequation 2996 9 \begin{equation} \label{eq106} \psi \equiv \frac{175}{3\pi ^2}A_{\rm h} (x)x^{-1/2}\Gamma . \end{equation}(C.14j)The above relations constitute the analytical solution of the second part of the problem (see C99, Figs. 1–16).

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.