A&A 384, 1119-1123 (2002)
DOI: 10.1051/0004-6361:20011773

Critical Richardson numbers and gravity waves

V. M. Canuto

NASA, Goddard Institute for Space Studies, 2880 Broadway, New York, NY 10025, USA
Dept. of Applied Physics and Mathematics, Columbia University, New York, NY 10027, USA

Received 23 July 2001 / Accepted 7 December 2001

Abstract
In this paper we present two new results. The first concerns the proper identification of the critical Richardson number Ri(cr) above which there is no longer turbulent mixing. Thus far, all studies have assumed that:

$\displaystyle % Ri {\rm (cr)} = Ri^\ell {\rm (cr)} = 1/4.$     (1)

However, since $Ri^\ell$(cr) determines the upper limit of a laminar regime (superscript $\ell$), it has little relevance to stars where the problem is not to determine the end point of a laminar regime but the endpoint of turbulence. We show that the latter is characterized by $Ri^{\rm t} {\rm (cr)}$, where t stands for turbulence, and has a value four times larger than (1):
$\displaystyle % Ri {\rm (cr)} = Ri^{\rm t} {\rm (cr)} \approx 4 Ri^\ell {\rm (cr)}\approx 1.$     (2)

We also show that use of (2) instead of (1) changes the conclusions of recent studies. Inclusion of radiative losses (characterized by the Peclet number Pe) which weaken stable stratification and help turbulence, further changes (2) to (r stands for radiative):
$\displaystyle % Ri {\rm (cr)} = Ri^{\rm r} {\rm (cr)} \sim (1+Pe)Pe^{-1}Ri^{\rm t} {\rm (cr)}$     (3)

which, for Pe<1, allows turbulence to survive far longer than (2). Finally, turbulent convection generates gravity waves that propagate into the radiative region and act as an additional source of energy. This further changes Eq. (3) to (gw stands for gravity waves):
$\displaystyle % Ri {\rm (cr)} = Ri^{\rm gw} {\rm (cr)} =Ri^{\rm r} (1+\eta_{\rm gw})$     (4)

where $\eta_{\rm gw}>1$. In conclusion, the successive inclusion of relevant physical processes leads to a chain of increasing values of $Ri {\rm (cr)}$:
$\displaystyle % Ri {\rm (cr)} = Ri^{\ell} {\rm (cr)} \rightarrow Ri^{\rm t} {\rm (cr)}\rightarrow Ri^{\rm r} {\rm (cr)}\rightarrow Ri^{\rm gw} {\rm (cr)}.$     (5)

The second result concerns the dependence of the diffusivity D on $\Omega $. We show that the commonly used expression
$\displaystyle % D_{\chi^{-1}}\sim (\Omega/N)^2$     (6)

is not correct for the regime Pe<1 that characterizes a stably stratified regime. The proper $\Omega $-dependence is:
$\displaystyle % D_{\chi^{-1}}\sim (\Omega/N)^4.$     (7)

Key words: Sun: general - Sun: interior - convection - turbulence


1 The problem

The determination of turbulent diffusivities for a stably stratified, radiative regime in the presence of sources like shear (differential rotation, meridional currents), encompasses several aspects that need to be analyzed. First, since a diffusivity D is of the form:
$\displaystyle % D\approx w\Lambda$     (1a)

where w is an rms velocity and $\Lambda$ is a mixing length, one needs to specify two independent variables, w and $\Lambda$. Using Kolmogorov law, $K\sim (\varepsilon \Lambda)^{2/3}$, where K is the turbulent kinetic energy and $\varepsilon$ is the rate of dissipation of K, Eq. (1a) can be written as:
$\displaystyle % D\approx K^2 \varepsilon^{-1}$     (1b)

which is physically much clearer than (1a). In fact, K has its maximum at the largest scales whereas $\varepsilon$ (being a dissipation) peaks at the smallest scales. The specification of K and $\varepsilon$ is therefore equivalent to specifying both large and small scales. Turbulence modeling provides the two dynamic equations satisfied by Kand $\varepsilon$ (Canuto 1998; Kupka 2001).

The second problem is that a large variety of independent data from laboratory to large eddy simulations (Webster 1964; Istweire & Helland 1989; Schumann & Gerz 1995; Wang et al. 1996; Kosovic & Curry 2000; Strang & Fernando 2001) have shown that stratification has different effects on heat and momentum transport but Eq. (1b) is unable to account for that fact. On the other hand, turbulence models yield the following result:

$\displaystyle % D_\alpha = 2S_\alpha K^2 \varepsilon^{-1}$     (1c)

where the $S_\alpha$ ( $\alpha = m,h$, momentum and heat) are dimensionless structure functions that are different from momentum and heat. The $S_\alpha$ are functions of:
$\displaystyle % S_\alpha = S_\alpha (Pe,\,Ri)$     (1d)

where Pe is the Peclet number that characterizes radiative losses ($\chi$ in cm2 s-1 is the radiative diffusivity):
$\displaystyle % Pe = \frac{4\pi^2}{125}\frac{K^2}{\varepsilon\chi}$     (1e)

Ri is the Richardson number:
$\displaystyle % Ri = N^2\Sigma^{-2}$     (2a)

characterizing the ratio between sinks (stable stratification N2 and sources (shear instabilities, $\Sigma^{2}$). Specifically, $\Sigma$ is the mean rate of strain ( $\bar{\vec u}$ is the mean velocity field)
$\displaystyle % \Sigma^{2} = 2\Sigma_{ij}\Sigma_{ij},\quad 2\Sigma_{ij} = \bar{u}_{i,j} + \bar{u}_{j,i}$     (2b)

and N is the Brunt-Vaisala frequency:
$\displaystyle % N^{2} = g\alpha \left[\frac{\partial T}{\partial z} - \left(\fr... ...rtial z}\right)_{\rm ad}\right] = gH_{\rm p}^{-1} (\nabla_{\rm ad} - \nabla)>0.$     (2c)

Third, there must exist a critical value Ri(cr) of the Richardson number at which turbulent mixing ceases since stable stratification has become too strong. Operationally, one defines $Ri^{\rm t}$(cr) as the value at which the turbulent kinetic energy vanishes:
$\displaystyle % K\left[Ri^{\rm t} {\rm (cr)}\right] = 0$     (2d)

and Eq. (1c) implies that the diffusivities also vanish at $Ri^{\rm t}$(cr). Given the above procedure and a turbulence model that provides the functions:
$\displaystyle % K,\,\varepsilon,\,S_\alpha$     (2e)

the value of $Ri^{\rm t}$(cr) is given by solving (2d) and need not be imposed. When no turbulence model is available and the diffusivities are constructed using heuristic arguments, the lack of a dynamic equation for K deprives such models of Eq. (2d) and $Ri^{\rm t}$(cr) must be guessed, a procedure that has been used thus far, as we discuss next.

2 The critical Richardson number

Lacking a turbulence model, in all astrophysical literature thus far (Zahn 1974; Pinsonneault et al. 1989; Meader 1995; Meader & Meynet 1996; Garcia-Lopez & Spruit 1991; Talon & Zahn 1997; Schatzman et al. 2000), it has been customary to identify:

$\displaystyle % Ri{\rm (cr)} = Ri^\ell {\rm (cr)}$     (3a)

where
$\displaystyle % Ri^\ell {\rm (cr)} = 1/4.$     (3b)

Equation (3b) was derived by Miles (1961) and Howard (1961) using linear stability analysis and represents the value at which a laminar flow (thus the superscript $\ell$) becomes unstable. In what follows, several arguments are presented to show that, given the specific task of identifying when, in a stably stratified region of a star, turbulent mixing ceases to exist, Eqs. (3a), (3b) are not the correct answer.

First, Eq. (3b) is the result of a "road from laminarity'', an approach that can only tell us when a linear regime breaks down and non-linearities becom important. When that occurs, one first enters a weakly non-linear regime and then finally in a turbulent regime where the non-linearities dominate. We refer the reader to an illustrative depiction of these successive processes (Woods 1969) which highlights the different physical regimes of laminar, weakly non-linear and strongly non-linear regimes. Given a stable laminar sheet of thickness h, Kelvin-Helmholtz instabilities gradually erode and entrain fluid parcels above and below h. The process leads to an incrase of h which ceases when the thickness has become about four times the original value h. Woods (1969) concluded that "since the final thickness is nearly four times the original value, the final Richardson number is also four times the value prior to the instability''. The turbulent Richardson number $Ri^{\rm t}$ is thus computed to be:

$\displaystyle % Ri^{\rm t} {\rm (cr)} \approx 4 Ri^\ell {\rm (cr)} \approx 1.$     (3c)

Equation (3c) has also been derived using other arguments. Once the regime described by (3b) has been crossed, non-linearities become important and one must carry out a stability analysis including them. Abarbanel et al. (1984) proved that in such a case the sufficient and necessary condition for instability is:
$\displaystyle % Ri\leq Ri^{\rm t} \approx 1$     (3d)

which agrees with the heuristic argument leading to (3c).

A more physical argument relies on energy considerations. A turbulent state can be maintained only as long as the sources outweigh the sinks. Using the equation for the turbulent kinetic energy (Canuto 1998), this criterion becomes:

$\displaystyle % R_{\rm f} <1$     (4a)

where $R_{\rm f}$ is the "flux Richardson number'' defined as (see Eqs. (9b)-(9d)):
$\displaystyle % R_{\rm f} = \frac{{\rm heat\ flux}}{{\rm momentum\ flux}}= \frac{D_{\rm h}}{D_{\rm m}} Ri <1.$     (4b)

Thus, Eq. (4a) becomes:
$\displaystyle % Ri <Ri^{\rm t} {\rm (cr)} \equiv \frac{D_{\rm m}}{D_{\rm h}} = \frac{S_{\rm m}}{S_{\rm h}} \cdot$     (4c)

Martin (1985) showed that the computed oceanic upper mixed layer (which is strongly stable and stirred by shear) could match observations only if (3d) was satisfied. Large eddy simulation of a stable planetary boundary layer by Kosovic & Curry (2000) also conclude that mixing exists up to $Ri^{\rm t}\approx 1$. A number of other cases leading to (3d) are discussed in Strang & Fernando (2001).

From the turbulence modeling viewpoint, turbulence models (Canuto 1998) provide the dependence of the turbulent kinetic energy K on Ri. Since the stronger the degree of stability, the harder it is for turbulence to survive, K is a decreasing function of Ri and at some $Ri^{\rm t} {\rm (cr)}$, Eq. (2d) is satisfied and the value of $Ri^{\rm t} {\rm (cr)}$ so identified. The solution of Eq. (2d) was obtained in (Canuto 1998, Fig. 6a) and the result agrees with (3c):

$\displaystyle % Ri^{\rm t} {\rm (cr)} \approx 1.$     (4d)

In conclusion, $Ri^{\ell} {\rm (cr)}$ and $Ri^{\ell} {\rm (cr)}$ represent two different physical situations. In laboratory experiments one follows the development of an instability and witnesses the transition from a laminar to a weakly non-linear state, a bottom-up approach that confirms $Ri^{\ell} {\rm (cr)}$ as the transition value. However, in stars one cannot follow such a development and one must therefore adopt a conceptually different approach, a top-down approach. One is presented with a situation in which there is turbulent mixing and the relevant question then changes to: what is maximum value of Ri for which such mixing can be sustained? Since turbulent mixing is due to strong non-linearities, the answer is $Ri^{\rm t} {\rm (cr)} \approx 4 Ri^{\ell} {\rm (cr)}$. Use of $Ri^{\ell} {\rm (cr)}$ to characterize the boundary between turbulent mixing and no mixing is incorrect since $Ri^{\ell} {\rm (cr)}$ refers to a physical situation quite different from the one most astrophysical studies intend to describe.

In conclusion, one should forgo the bottom-up approach leading to $Ri^{\ell} {\rm (cr)}$ and adopt the top-down approach that leads to $Ri^{\rm t} {\rm (cr)}$.

3 Two consequences

At first sight, the change from 1/4 to 1 may seem marginal but it is not so, as we now discuss. Garcia-Lopez & Spruit (1991) showed that for gravity waves to be a source of mixing compatible with the data on Li7 depletion, their result had to be boosted by an arbitrary factor of ${\approx}15$ which meant that gravity waves mixing is inefficient. However, had the authors used $Ri^{\rm t} {\rm (cr)}$ instead of $Ri^{\ell} {\rm (cr)}$, the boosting factor would have been only ${\approx}3$, which is well within the uncertainties of the problem. Their conclusion would have been that gravity waves are an efficient source of mixing. More recently, Schatzman et al. (2000, SZM) studied whether mixing may be due to shear instabilities. They concluded that for that to occur (their Eq. (8)):

$\displaystyle % \xi > 40/Ri{\rm (cr)},\quad \xi \equiv (\chi/\nu) Ri^{-1}$     (5a)

where $\nu$ and $\chi$ are the molecular viscosity and radiative conductivity. SZM used:
$\displaystyle % Ri{\rm (cr)} = Ri^{\ell} {\rm (cr)} = 1/4$     (5b)

and showed (their Fig. 2) that $\xi$ almost never satisfies (5a). Had SZM used:
$\displaystyle % Ri{\rm (cr)} = Ri^{\rm t} {\rm (cr)} = 1$     (5c)

the first relation in (5a) would have been satisfied and SZM would have concluded that shear instabilities do set in. Whether shear instabilities are sufficiently strong to explain the diffusiviy required to explain the depletion of Li7 is a different matter since the quantification of how much mixing is produced by an instability requires a model much more complex than linear analysis, that is, a turbulence model.

4 Peclet number dependence

As we stressed earlier, the solution of the Reynolds stress equations provided by a turbulence models (Canuto 1998) yields not only the form (1c) but the structure functions $S_\alpha$ as well. Equation (4c) can therefore be written as (r stands for radiative losses):

$\displaystyle % Ri<Ri^{\rm r}{\rm (cr)}$     (6a)


$\displaystyle % Ri^{\rm r}{\rm (cr)}= \frac{S_{\rm m} (Ri,Pe)}{S_{\rm h} (Ri,Pe)}\cdot$     (6b)

The Pe-dependence of the functions $S_\alpha$ is such that:
$\displaystyle % Pe >1{:}\, S_{\rm m} (Ri,Pe)\approx S_{\rm m} (Ri),\ \ S_{\rm h}(Ri,Pe)\approx S_{\rm h}(Ri).$     (7a)

Both $S_{\rm m,h}$ are independent of Pe. In the opposite case:
$\displaystyle % Pe <1{:}\, S_{\rm m} (Ri,Pe){\approx} S_{\rm m} (Ri),\ \ S_{\rm h}(Ri,Pe){\approx }Pe S_{\rm h}(Ri)$     (7b)

$S_{\rm m} (Ri,Pe)$ is still Pe-independent while $S_{\rm h}(Ri,Pe)$ depends linearly on Pe. Combining (7a,b) into a single expression and substituting into (6b), we obtain:
$\displaystyle % Ri^{\rm r}{\rm (cr)}= \frac{1 + Pe}{Pe}\frac{S_{\rm m} (Ri)}{S_{\rm h}(Ri)} = \frac{1+Pe}{Pe} Ri^{\rm t}{\rm (cr)}$     (8a)

which, in the limit of small Pe, yields:
$\displaystyle % Ri^{\rm r}{\rm (cr)}\approx Pe^{-1}Ri^{\rm t}{\rm (cr)}.$     (8b)

Thus, $Ri^{\rm r}{\rm (cr)}$ can be much larger than $Ri^{\rm t} {\rm (cr)}$. Using (8a), we rewrite (6a) as:
$\displaystyle % Pe <1\qquad Ri < \frac{1}{Pe}\frac{S_{\rm m} (Ri)}{S_{\rm h}(Ri)}\cdot$     (8c)

5 Gravity waves

In principle, once we add to the previous formalism the dynamic equations for K and $\varepsilon$, the model is complete. The results can be found in Canuto (1998). However, such an approach, while complete from the point of view of turbulence modeling, does not take into account the presence of gravity waves which are generated in the convective zone and propagate into the radiative zone with a flux $\Pi$(gw). Such waves act a source of turbulence (in addition to that of shear e.g., differential rotation and meridional currents). In the solar case, Kumar et al. (1999) have computed the flux of gravity waves to be:

$\displaystyle % \Pi{\rm (gw)} = 0.1\% (L/M)_{\rm sun} = 2\times 10^{-3}~{\rm cm^2~s^{-3}}.$     (9a)

The question then arises: how can one include $\Pi$ into a turbulence model? Here, we suggest a first approach. Consider the stationary limit of the time dependent equation for the turbulent kinetic energy (Canuto 1998, Eq. (25a)):
$\displaystyle % {\rm Production = Dissipation}$     (9b)

where:
$\displaystyle % {\rm Production} = - \tau_{ij}\Sigma_{ij} + g\alpha {\overline{w\theta}} + \Pi.$     (9c)

The first term represents the source due to shear (e.g., differential rotation, meridional currents), the second term is the sink due to stable stratification and the third term is the contribution of gravity waves. Using the following forms for the Reynolds stresses $\tau_{ij}$ and the heat flux $\overline{w\theta}$ (Canuto 1998, Eqs. (39a) and (43a)):
$\displaystyle % \tau_{ij}= - 2D_{\rm m} \Sigma_{ij},\quad \overline{w\theta} = ... ...}{\partial z} - \left(\frac{\partial T}{\partial z}\right)_{\rm ad}\right]\cdot$     (9d)

Equations (9b), (9c) become:
$\displaystyle % D_{\rm m} \Sigma^2 - D_{\rm h}N^2 + \Pi = \varepsilon$     (10a)

which can be written as:
$\displaystyle % \varepsilon = (D_{\rm m} \Sigma^2 + \Pi)(1 - R_{\rm f}^{\rm gw})$     (10b)

where the gravity wave-dependent flux Richardson number is defined as:
$\displaystyle % R_{\rm f}^{\rm gw}\equiv \frac{D_{\rm h}}{D_{\rm m}}\frac{Ri}{1... ...m}\Sigma^2} = \frac{S_{\rm h}}{S_{\rm m}}\frac{Ri}{1 + \Pi/D_{\rm m}\Sigma^2}<1$     (10c)

which replaces (4b). Using (6b), we obtain:
$\displaystyle % Ri < Ri^{\rm gw}{\rm (cr)}$     (10d)

where:
$\displaystyle % Ri^{\rm gw}{\rm (cr)} = Ri^{\rm r}{\rm (cr)} (1 + \Pi/D_{\rm m}\Sigma^2)>Ri^{\rm r}{\rm (cr)}.$     (10e)

Thus, in the presence of gravity waves, Rigw(cr) is larger than Rir(cr).

6 $\Omega $ dependence of the scalar diffusibity

The expression for the scalar turbulent diffusivity D employed thus far in the literature does not differentiale between momentum and heat. Its specific form has been taken to be:

$\displaystyle % Pe < 1\qquad D{\chi^{-1}}\sim \Omega^2 N^{-2}$     (11a)

where and $\Omega^2$ stands for $\Sigma^{2}$. Equation (11a) has been widely used (Baglin et al. 1985; Maeder 1995; Maeder & Meynet 1996; Talon & Zahn 1997; Charbonnel & Talon 1999).

Let us repeat the standard derivation of (11a) keeping the $S_\alpha$ structure functions. Starting with (1c) and using (1e), one has $(a\equiv 125/2\pi^2)$:

$\displaystyle % D_\alpha \chi^{-1}= aPe\,S_\alpha (Ri,Pe).$     (11b)

For momentum and heat diffusivities, we have:
$\displaystyle % D_{\rm m} \chi^{-1}$ = $\displaystyle aPe\,S_{\rm m} (Ri,Pe)$ (11c)
$\displaystyle D_{\rm h} \chi^{-1}$ = $\displaystyle aPe\,S_{\rm h} (Ri,Pe).$ (11d)

Because of (7a), (7b), we can further write for any Pe:
$\displaystyle % D_{\rm m} \chi^{-1}\sim aPe\,S_{\rm m} (Ri)$     (11e)
$\displaystyle D_{\rm h} \chi^{-1}= aPe^2 (1+Pe)^{-1}\,S_{\rm h} (Ri).$     (11f)

The different Pe-dependence of $D_{\rm m}$ and $D_{\rm h}$ is a result that was first derived by Howells (1960) using heuristic arguments and later proven rigorously with RNG (Renormalization group; Canuto & Dubovikov 1996, Eqs. (C16-18)). Next, we employ (8c) to eliminate Pe and obtain:
$\displaystyle % Pe < 1\qquad D_{\rm m} \chi^{-1}\sim S(Ri)(\Omega/N)^2$     (12a)
$\displaystyle D_{\rm h} \chi^{-1}\sim S(Ri)(\Omega/N)^4$     (12b)

where the function S(Ri)
$\displaystyle % S(Ri)\equiv S^2_{\rm m}(Ri)/S_{\rm h}(Ri)$     (12c)

is plotted in Fig. 1. When compared with (11a), the new expressions (12) exhibit several new features:
1)
heat and momentum diffusivity are different,
2)
only the momentum diffusivity scales like $\Omega^2/N^2$,
3)
heat diffusivity scales like $(\Omega^2/N^2)^2$,
4)
the absence of the $S_\alpha$ functions in phenomenological models means that they lack the function S(Ri) which ensures that the diffusivities vanish at some point (Fig. 1),
5)
the value at which $D_{\rm m,h}$ vanish is represented by the variable S(Ri), see Fig. 1. The point at which the diffusivities vanish corresponds to the value of Ri at which the turbulent kinetic energy vanishes.
In summary, the above derivation shows the major differences brought about by the structure fonctions $S_\alpha$ vis à vis the standard model where they are absent.
  \begin{figure} \par\includegraphics[width=6.3cm,clip]{fig1.eps}\end{figure} Figure 1: The function S(Ri) defined in Eq. (12c) is plotted vs. the Richardson number Ri. As one can observe, turbulence exists way past the Ri= 1/4 value. It dies only at $Ri \sim 1$.
Open with DEXTER

7 Conclusions

Several arguments have been presented to show that in a stably stratified fluid, shear dominated turbulence persists above the laminar value $Ri^\ell {\rm (cr)} = 1/4$. Regrettably, the lack of a turbulence model has prompted the use of the laminar value $Ri^{\ell} {\rm (cr)}$ in spite of the fact that it does not represent the physical situation encountered in stellar interiors.

The bottom-up approach starting from linear stability and leading to $Ri^\ell {\rm (cr)} = 1/4$ must be substituted with a top-down approach that gives the range of values of Ri for which turbulent mixing operates. The new $Ri^{\rm t} {\rm (cr)}$ is defined as the value at which turbulence ceases to exist since the turbulent kinetic energy vanishes. Several independent arguments are presented to show that:

$\displaystyle % Ri^{\rm t} {\rm (cr)}\approx 1$     (13a)

a value that would have changed the conclusions of Schatzman et al. (2000) and considerably strengthen the conclusions of Garcia-Lopez & Spruit (1991) concerning gravity waves as a source of mixing.

Radiative losses further increase the critical $Ri^{\rm t} {\rm (cr)}$ to $Ri^{\rm r}{\rm (cr)}$ which for Pe <1 is ${\approx} Pe^{-1}$ times larger than (13a). Gravity waves act as a source of mixing leading to a further boost in the value of $Ri {\rm (cr)}$. Through the successive inclusion of physical processes we have shown the sequence:

$\displaystyle % Ri^\ell {\rm (cr)}\rightarrow Ri^{\rm t} {\rm (cr)}\rightarrow Ri^{\rm r} {\rm (cr)}\rightarrow Ri^{\rm gw} {\rm (cr)}$     (13b)

with
$\displaystyle % Ri^\ell {\rm (cr)}< Ri^{\rm t} {\rm (cr)}< Ri^{\rm r} {\rm (cr)}< Ri^{\rm gw} {\rm (cr)}.$     (13c)

In conclusion, the traditional $Ri^\ell {(cr)} = 1/4$ must be replaced by Ri gw (cr).

In many ways, the new formalism resolves the dichotomy that has existed thus far between shear-dominated and gravity wave-dominated scenarios. Both processes have been integrated into a single formalism. Numerical results leading to numerical values of the scalar diffusivity D are presented in an upcoming paper (Canuto & Minotti 2001).

Acknowledgements
The author would like to thank an anonymous referee for inquisitive questions that helped improve the original manuscript.

References

 
Copyright ESO 2002