# Turbulent Transport by Diffusive Stratified Shear Flows: From Local to Global Models. I. Numerical Simulations of a Stratified Plane Couette Flow

Pascale Garaud, Damien Gagnier, Jan Verhoeven

► To cite this version:

Pascale Garaud, Damien Gagnier, Jan Verhoeven. Turbulent Transport by Diffusive Stratified Shear Flows: From Local to Global Models. I. Numerical Simulations of a Stratified Plane Couette Flow. The Astrophysical Journal, 2017, 837, 10.3847/1538-4357/837/2/133. insu-03677072

HAL Id: insu-03677072

https://insu.hal.science/insu-03677072v1

Submitted on 23 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Turbulent Transport by Diffusive Stratified Shear Flows: From Local to Global Models.
## I. Numerical Simulations of a Stratified Plane Couette Flow

Pascale Garaud¹, Damien Gagnier², and Jan Verhoeven³

¹Department of Applied Mathematics and Statistics, Baskin School of Engineering, University of California at Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

²Institut de Recherche en Astrophysique et Planétologie (IRAP), 14, avenue Edouard Belin, F-31400 Toulouse, France

³Department of Earth and Planetary Sciences, University of California at Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA
Received 2016 October 12; revised 2016 December 10; accepted 2016 December 12; published 2017 March 10

## Abstract

Shear-induced turbulence could play a significant role in mixing momentum and chemical species in stellar radiation zones, as discussed by Zahn. In this paper we analyze the results of direct numerical simulations of stratified plane Couette flows, in the limit of rapid thermal diffusion, to measure the turbulent viscosity and the turbulent diffusivity of a passive tracer as a function of the local shear and the local stratification. We find that the stability criterion proposed by Zahn, namely that the product of the gradient Richardson number and the Prandtl number must be smaller than a critical values (JPr)$_c$ for instability, adequately accounts for the transition to turbulence in the flow, with (JPr)$_c$ ≈ 0.007. This result recovers and confirms the prior findings of Prat et al. Zahn's model for the turbulent diffusivity and viscosity, namely that the mixing coefficient should be proportional to the ratio of the thermal diffusivity to the gradient Richardson number, does not satisfactorily match our numerical data. It fails (as expected) in the limit of large stratification where the Richardson number exceeds the aforementioned threshold for instability, but it also fails in the limit of low stratification where the turbulent eddy scale becomes limited by the computational domain size. We propose a revised model for turbulent mixing by diffusive stratified shear instabilities that properly accounts for both limits, fits our data satisfactorily, and recovers Zahn's model in the limit of large Reynolds numbers.

**Key words:** hydrodynamics – instabilities – stars: general – turbulence

## 1. Introduction

Thanks to recent advances in supercomputing it is now possible to run numerical experiments designed to quantify the rates of turbulent mixing of selected quantities (heat, angular momentum, composition, etc.) in the presence of various hydrodynamic and magnetohydrodynamic instabilities, in parameter regimes appropriate for astrophysical plasmas that cannot be achieved in more traditional laboratory experiments. This raises the engaging prospect of finally being able to constrain models of non-canonical mixing in stars through first-principles theory and experiments rather than through observations alone.

A strong potential candidate for driving vertical mixing in stellar radiation zones is the shear instability, which has gained significant popularity since the discovery of relatively strong radial shear layers via helio- and asteroseismology, in the Sun (Christensen-Dalsgaard & Schou 1988; Brown et al. 1989) and in Red Giant Branch stars (Deheuvels et al. 2012, 2014). However, the radiation zones of these stars are usually too strongly stratified to be shear-unstable, at least if thermal diffusion is ignored. Indeed, the standard criterion for instability of inviscid, non-diffusive stratified shear flows to infinitesimal perturbations is the Richardson criterion (formally proved by Howard 1961; Miles 1961), which states that the local gradient Richardson number

$$J = \frac{N^2}{S^2}, \qquad (1)$$

where N is the Brunt-Väisälä frequency and S is the local shearing rate, must drop below 1/4 somewhere in the flow for instability. It is sometimes argued that finite-amplitude non-

diffusive instabilities may exist for somewhat larger values of J, but energetic arguments (Richardson 1920) suggest that J cannot be much larger than 1 for sustained turbulent motions to exist.

In thermally stratified systems, Townsend (1958) and Dudis (1974) noted, however, that thermal radiation and/or thermal diffusion let the perturbed fluid adjust thermally to its surroundings, thereby reducing its stabilizing buoyancy excess or deficit compared with the background as it moves up and down. In an optically thick fluid, the ratio of the thermal diffusion timescale to the shearing timescale for a shear layer of vertical lengthscale L is measured by the global Péclet number,

$$Pe_L = \frac{SL^2}{\kappa_T}, \qquad (2)$$

where $\kappa_T$ is the thermal diffusivity. Linear stability analyses have shown that having $Pe_L \ll 1$ can destabilize stratified shear layers that would otherwise be stable in the non-diffusive case (Jones 1977; Lignières et al. 1999; Garaud et al. 2015), and typically raises the critical Richardson number for instability by a factor proportional to $Pe_L^{-1}$. Note that unfortunately there is no simple formal criterion for linear instability equivalent to the Miles-Howard theorem when thermal diffusion is accounted for, even in the inviscid limit. Also note that, in the optically thick limit, there is no geophysical analog of turbulent and strongly diffusive stratified shear flows. This is because in geophysical fluids the Prandtl number $Pr = \nu/\kappa_T$ (where $\nu$ is the kinematic viscosity) is always of order one or larger, and since the global Reynolds number $Re_L \equiv SL^2/\nu$, we have $Re_L = Pe_L/Pr < Pe_L \ll 1$. In other words, low Péclet number
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of the global vs. local properties of stratified shear flows. The shading represents the stratification, from denser fluid (blue) to lighter fluid (red). On the global scale, the shear may not be constant but will have a well-defined amplitude *U* and lengthscale *L*, defining a global Péclet number Pe<sub>*l*</sub> = *UL*/κ<sub>*T*</sub>, which is likely much larger than one. On small enough scales *l*, however, one can define an eddy Péclet number Pe<sub>*l*</sub> = *Sl*<sup>2</sup>/κ<sub>*T*</sub> that is much smaller than 1. On these small scales, the local shearing rate *S* is approximately constant. Zahn's models assume that the turbulence develops on these small vertical scales, so both stability criterion and mixing can be treated in a local fashion.</figcaption></figure>

flows in geophysics always have a low Reynolds number as well, implying that their dynamics must be viscously dominated. In stellar astrophysics, however, Pr ≪ 1 so it is possible to have, at the same time, Pe<sub>*l*</sub> ≪ 1 and Re<sub>*l*</sub> ≫ 1. This limit is unique to astrophysical applications, so none of the extensive literature currently available on geophysical stratified shear flows applies.

In any case, linear theory is of limited use, since it neither addresses the possibility of destabilization by finite-amplitude perturbations, nor the question of the saturation of the instability, which is necessary to understand the properties of the turbulent flows that later develop. In an attempt to go beyond linear theory, Zahn therefore proposed a now widely used two-part model for shear instabilities in stars that involves a revised criterion for instability (Zahn 1974) and a turbulent mixing prescription (Zahn 1992). We describe both in turn.

1.1. *Stability Criterion*

To address the question of stability, Zahn (1974) noted that, allowing for finite-amplitude instabilities, structures of height *l* much smaller than the shear scaleheight *L* can in principle be destabilized provided their eddy Péclet number Pe<sub>*l*</sub> ≡ *Sl*<sup>2</sup>/κ<sub>*T*</sub> is much smaller than 1 (even if Pe<sub>*l*</sub> ≫ 1); see Figure 1. By analogy with the results of Townsend (1958) and Dudis (1974), he proposed that instability on the scale *l* may occur, in the inviscid limit, as long as

$$
JPe_l < (JPe)_c, \tag{3}
$$

where (JPe)$_c$ is a constant of order one.

Accounting for the effect of viscosity limits the range of
scales that can be destabilized this way, requiring in addition
that

$$
Re_l \equiv \frac{Sl^2}{\nu} > Re_c. \quad (4)
$$

The constant Re<sub>c</sub> is a critical Reynolds number for finite-amplitude instability, which he picks to be of order 10<sup>3</sup>. Combining these two criteria leads to Zahn's local criterion for diffusive stratified shear instabilities (Zahn 1974), namely that

a flow is unstable provided

$$
J\mathrm{Pr} = J \frac{\mathrm{Pe}_l}{\mathrm{Re}_l} < \frac{(J\mathrm{Pe})_c}{\mathrm{Re}_c} \equiv (J\mathrm{Pr})_c. \quad (5)
$$

This defines the critical Richardson–Prandtl number for shear instabilities (JPr)$_c$, which, according to the arguments given above (Zahn 1974), should be of order $10^{-3}$. Note that, by contrast with (3), (5) is now independent of the turbulent eddy scale, which is quite convenient since the latter is not known a priori.

While Zahn's criterion is based on somewhat heuristic arguments, two important pieces of formal theory support it. First, as established by Lignières (1999), the dynamics of strongly diffusive stratified shear instabilities, that is, in the limit of asymptotically low Péclet numbers, can only depend on two non-dimensional parameters (aside from possible geometric effects), namely the global Reynolds number Re<sub>L</sub> and the product JPe<sub>L</sub>. Hence any stability criterion, should it exist, must be expressed in terms of Re<sub>L</sub> and JPe<sub>L</sub>, rather than Re<sub>L</sub>, J and Pe<sub>L</sub> separately. It is easy to verify that it is indeed true in Zahn's criterion since JPr = JPe<sub>L</sub>/Re<sub>L</sub>.

Second, in this low Péclet number limit (Lignières 1999) it is
possible to obtain formal results on the energy stability of
stratified shear flows. Energy stability is arguably a more useful
concept than linear stability because it implies that any initial
perturbation must decay, regardless of their amplitude or form,
while linear stability only implies that infinitesimal perturba-
tions must decay. As a result, establishing the energy stability
of a shear layer guarantees that there can be no shear-induced
turbulent mixing (unless other physics are added that could
further destabilize it). For instance, Bischoff (2013) showed
that a low Péclet number linear shear in a linearly stratified
background (i.e., where *S* and *N* are both constant) is energy
stable provided *JPr* > 1/4. A similar energy stability criterion
was obtained by Garaud et al. (2015) for any periodic shear
flow, albeit with a constant that is different from 1/4 and that
depends on the shape of the domain considered. One may
therefore postulate that a reasonable general energy stability
criterion for diffusive stratified shear layers in stellar interiors is

$$
J\mathrm{Pr} > (J\mathrm{Pr})_e, \tag{6}
$$

where (JPr)<sub>e</sub> is a constant of order unity that may depend on the
aspect ratio of the domain. The scaling of this criterion with *J*
and *Pr* is identical to the one proposed by Zahn (1974),
suggesting that his physically motivated derivation could be a
plausible explanation for the more formal but significantly less
intuitive energy stability criterion of Bischoff (2013) and
Garaud et al. (2015).

Despite these supporting arguments, however, neither linear
theory nor energy stability theory can really address the
question of whether Zahn's criterion is indeed a good predictor
for the transition to turbulence in stratified shear flows. A
linearly stable shear flow may still be the seat of vigorous
turbulent mixing. Meanwhile, a shear layer lying outside of the
energy-stable region of parameter space does not necessarily
imply turbulence. What really matters in terms of stellar mixing
is neither energy stability nor linear instability, but instead the
more elusive threshold for sustained turbulence driven by
finite-amplitude instabilities. This threshold can only be
obtained using 3D direct numerical simulations (DNSs).
## 1.2. Mixing Model

Later, Zahn expanded his theory for stratified shear instabilities in stellar interiors, presenting a model for the induced turbulent diffusivity $D_{turb}$ (Zahn 1992). Starting from dimensional analysis, he first proposed that

$$D_{\text{turb}} = \beta S l_e^2, \qquad (7)$$

where $l_e$ is the lengthscale of the dominant turbulent eddies and $\beta$ is a constant of order unity. He then argued that $l_e$ can be estimated by assuming it is also the scale of turbulent eddies that are just marginally unstable according to criterion (3), or in other words, the largest lengthscale for which (3) holds:

$$J \frac{S l_e^2}{\kappa_T} = (J\text{Pe})_c, \qquad (8)$$

which implies

$$l_e = \sqrt{\frac{(J\text{Pe})_c \kappa_T}{JS}} = \sqrt{\frac{(J\text{Pe})_c}{J\text{Pe}_L}} L \equiv l_Z. \qquad (9)$$

We shall henceforth refer to this scale as the Zahn scale, and denote it as $l_Z$. Using $l_Z$ in $D_{turb}$, we arrive at

$$D_{\text{turb}} = \beta (J\text{Pe})_c \frac{\kappa_T}{J} \equiv C \frac{\kappa_T}{J} \quad \text{where} \quad C = \beta (J\text{Pe})_c \qquad (10)$$

is a universal constant.

Two points are worth noting. First, this model does not specify what quantity $D_{turb}$ actually mixes, and in this sense, could be applied both to the transport of chemical species or of momentum, albeit perhaps with a different constant $\beta$. Second, this model must clearly fail both in the limit of very large stratification (i.e., when the system is stabilized for $J\text{Pr} > (J\text{Pr})_c$ so $D_{turb}$ must drop to zero) and in the limit of very low stratification (where $J$ tends toward 0 and $D_{turb}$ would formally tend toward infinity, which is not physical). Whether the model holds between these two extreme limits, however, remains to be determined using 3D DNSs.

## 1.3. Numerical Simulations of Diffusive Stratified Shear Instabilities

Numerical simulations are therefore required to test the two parts of Zahn's model. The first three-dimensional DNSs of stratified shear instabilities in the diffusive limit (more specifically, in the limit of low Péclet number) were presented by Prat & Lignières (2013; see also Prat & Lignières 2014), followed by Garaud et al. (2015), Garaud & Kulenthirarajah (2016), and Prat et al. (2016). In all cases, they apply codes that use a horizontally periodic Cartesian domain, under the Boussinesq approximation. In nearly all cases (except in parts of the work of Prat & Lignières 2014), chemical species are mixed as passive tracers and stratification is due to temperature alone, which is also the limit we consider here. The background temperature gradient is assumed to be constant (in the low Péclet number limit), and therefore so is the buoyancy frequency $N$.

The works of Prat and his collaborators use two different model setups in which the imposed shear flow is linear everywhere in the domain, with a given constant shearing rate $S$. By doing so, they are testing Zahn's model as it was originally intended, namely as a local mixing model. To drive this linear shear flow, Prat & Lignières (2013, 2014) add a body force to

the momentum equation that takes the form of a relaxation term driving the horizontal average of the streamwise component of the velocity field toward the linear target flow $S_z$, where $z$ is the vertical coordinate, on a certain timescale $\tau$. Prat & Lignières (2014) argued that the results are independent of the exact value of $\tau$ provided the latter is small enough. More recently, Prat et al. (2016) revisited the homogeneously sheared model using the shearing sheet formalism instead. Crucially, they found that their results are consistent with those of Prat & Lignières (2013, 2014), hence validating the results of these studies, and showing that the manner in which the constant shear is forced does not influence the turbulent dynamics induced.

Prat et al. (2016) find that, in the low Péclet number limit, finite-amplitude instabilities exist up to about $J\text{Pr} \approx 0.007$, but initial perturbations decay for larger Richardson-Prandtl numbers. This constrains one of the unknown model constants, namely

$$(J\text{Pr})_c \approx 0.007, \qquad (11)$$

a value that is quite close to the one predicted by Zahn (1974), namely $(J\text{Pr})_c = O(10^{-3})$. They also measure the turbulent mixing coefficient $D_{turb}$ and use it to test the validity of the mixing prescription proposed in (10). They find that $D_{turb}J/\kappa_T$, which should be constant and equal to $C$, actually varies quite significantly with $J\text{Pr}$, and somewhat with Reynolds number $Re_L$. To be more specific, $D_{turb}J/\kappa_T$ first gradually increases with $J\text{Pr}$, reaches a maximum, then drops sharply to zero as $J\text{Pr} \to (J\text{Pr})_c$ (see Figure 6(a) in this paper for an equivalent figure made with our data). Meanwhile at low stratification $D_{turb}J/\kappa_T$ increases with $Re_L$, but becomes independent of it at larger stratification.

While the sharp drop in $D_{turb}J/\kappa_T$ as $J\text{Pr} \to (J\text{Pr})_c$ simply reflects the expected stabilization of the system against turbulence, Prat et al. (2016) have more difficulty explaining the decrease of $D_{turb}J/\kappa_T$ as $J\text{Pr}$ decreases. Looking at the turbulent eddies, however, they find that their vertical size increases as $J\text{Pr}$ decreases, and eventually approaches the domain size, at which point one should indeed expect Zahn's local model (10) to fail. They argue this also explains the observed dependence of $D_{turb}J/\kappa_T$ on the Reynolds number in the limit of low stratification. Indeed, since a larger Reynolds number corresponds to an intrinsically larger domain, it is possible to run simulations at lower values of $J\text{Pr}$ without being constrained by the domain size. Extrapolating their results to much larger Reynolds numbers (which are more appropriate of stellar conditions), Prat et al. (2016) then argue that $D_{turb}J/\kappa_T$ should in reality tend to a constant in the limit of very low stratification. They conclude by proposing the following prescription as a refinement of Zahn's model that fits their data in the limit of large Reynolds number:

$$\frac{D_{\text{turb}}J}{\kappa_T} \approx \alpha_0 + \alpha_1 J\text{Pr} - \alpha_2 (J\text{Pr})^2, \qquad (12)$$

where $\alpha_0 = 3.34 \times 10^{-2}$, $\alpha_1 = 18.8$ and $\alpha_2 = 2.86 \times 10^3$ (Prat et al. 2016).

In an attempt to move beyond homogeneous shear flows and investigate more complex shear profiles, Garaud and her collaborators (Garaud et al. 2015; Garaud & Kulenthirarajah 2016) chose to use a different model setup in which the shear is driven by a time-independent but sinusoidally varying body force $F(z) = F_0 \sin(kz)e_x$, where $k$ is a vertical wavenumber and $F_0$
is the amplitude of the force. By contrast with the work of Prat & Lignières (2013, 2014), and Prat et al. (2016), the shear is not specified in their simulations but is instead free to evolve in time either viscously or in response to the Reynolds stresses produced in the turbulent shear flow. As a result, once a quasi-stationary state has been achieved, the mean shear S varies with z and so does the gradient Richardson number J(z).

Measuring J(z) deep in the middle of the turbulent layer, Garaud & Kulenthirarajah (2016) found that turbulent solutions exist up to JPr ≃ 0.006, but not for significantly higher values of JPr. Furthermore, they found that the mixing coefficient (10) proposed by Zahn (1992) correctly models the transport of a passive scalar for a range of stratification, but fails either when it is too strong (due to the stabilization of the flow) or when it becomes too weak. In that case, they find that the size of the turbulent eddies becomes controlled by the lengthscale of the shear k⁻¹. All of these findings are entirely consistent with those of Prat et al. (2016).

A new coherent picture of diffusive stratified shear instabilities thus emerges from the combination of all these results. The stability criterion of Zahn (1974), given in (5), was verified to hold, and the constant (JPr)$_c$ was measured in several independent ways to be (JPr)$_c$ ≃ 0.007. The mixing prescription of Zahn (1992), given in (10), also appears to hold for intermediate values of the stratification. However, it fails in the limit of large stratification where the flow is stable, and in the limit of very low stratification where the turbulent eddy sizes become constrained either by the size of the domain (in the homogeneously sheared case) or by the shear lengthscale— in other words, when the local assumption at the heart of Zahn’s models fails.

A few loose ends remain to be investigated, however. First, while the limit of large stratification is reasonably well understood, the limit of low stratification remains to be studied in more detail. Is it possible, for instance, to extend Zahn’s mixing model (Zahn 1992) to account for it as well? Second, both Garaud & Kulenthirarajah (2016) and Prat et al. (2016) attempted to test Zahn’s mixing prescription more directly by measuring the dominant turbulent eddy scale $l_e$ and comparing $SI_e^2$ with $D_{\text{turb}}$ in the simulations. Both of these quantities are measured entirely independently (the first one being extracted from the power spectrum of the velocity fluctuations and the second one being extracted from the correlation between the velocity and passive tracer fields), and should be proportional to one another according to (7) for Zahn’s model to be correct. Curiously, Garaud & Kulenthirarajah (2016) and Prat et al. (2016) both found that $SI_e^2$ and $D_{\text{turb}}$ are not proportional to one another, even in the case of intermediate stratification where (10) roughly holds. This raises a strange conundrum (Garaud & Kulenthirarajah 2016): how can a model hold if its fundamental assumption does not?

To address these remaining questions and look at the problem from yet a new angle, we study DNSs of stratified plane Couette flows (i.e., stratified shear flows driven by plane-parallel no-slip plates moving in opposite directions) in the low Péclet number limit. In these numerical experiments, as we shall demonstrate, the shear flow between the plates is free to adjust itself (as in the body-forced case of Garaud and collaborators), but remains close to being linear in the bulk of the domain (as in the work of Prat and collaborators). This setup offers the advantage of being able to compare our results more easily with those of Prat et al. (2016) than with the

sinusoidally forced case. It is also motivated physically in the astrophysical context by the fact that many shear layers are forced by surface stresses rather than body forces. For instance, magnetic braking only affects the outer convective layers of a star, which then communicate the spin-down torque to the underlying radiative layer through interfacial stresses at the base of the convection zone. As another example, mass accretion onto a star, whether from a disk or a companion, usually deposits a thin layer of rapidly rotating fluid on the surface, which can also be viewed as imposing a surface stress to the underlying regions.

In what follows, we therefore present the stratified plane Couette flow setup in Section 2, and analyze the simulations in Section 3. We then revisit Zahn’s two-part model (Zahn 1974, 1992) in Section 4, and propose a new physically motivated mixing prescription for turbulent mixing by stratified shear instabilities that fits our numerical data in all limits. We conclude in Section 5 by summarizing our results and discussing avenues for future work.

## 2. The Model

The plane Couette flow is perhaps the simplest possible model setup for boundary-forced shear flows. Two no-slip horizontal infinite parallel plates, separated by a distance L, move opposite of one another other with velocities ±1/2 $\Delta U e_x$, driving a shear flow between them. If this flow remains entirely laminar, then it takes the form $\bar{U}_C(z) = \frac{\Delta U}{L}ze_x$, with a constant shearing rate $S_C = \Delta U/L$, the subscript C standing for “Couette.” If we further assume that the top and bottom plates are held at constant temperatures $T_m + \Delta T/2$ and $T_m - \Delta T/2$, in the absence of temperature fluctuations, the background temperature gradient is constant and equal to $T_{0z} = \Delta T/L$. If in addition $T_{0z} > T_z^{\text{ad}}$ where $T_z^{\text{ad}}$ is the adiabatic temperature gradient, this system is stably stratified. The dimensional equations and boundary conditions describing the evolution of the velocity field $\mathbf{u}$ and temperature fluctuations $T$ away from the background temperature profile $T_0(z) = T_m + T_{0z}z$, in the Boussinesq limit (Spiegel & Veronis 1960), are

$$ 
\begin{aligned}
\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} &= -\frac{1}{\rho_m} \nabla p + \alpha g T e_z + \nu \nabla^2 \mathbf{u}, \\
\nabla \cdot \mathbf{u} &= 0, \\
\frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T + w (T_{0z} - T_z^{\text{ad}}) &= \kappa_T \nabla^2 T,
\end{aligned}
\quad (13) $$

and

$$ 
\begin{aligned}
&v = w = T = 0 \text{ at } z = \pm L/2, \\
&u = \pm \Delta U/2 \text{ at } z = \pm L/2,
\end{aligned}
\quad (14) $$

together with periodic boundary conditions in the horizontal direction. In these equations, $\mathbf{u} = (u, v, w)$ is the full 3D velocity field, $p$ is the pressure perturbation away from hydrostatic equilibrium, $\rho_m$ is the mean density of the fluid, $\alpha = -\rho_m^{-1}(\partial\rho/\partial T)_p$ is the coefficient of thermal expansion, and $g$ is gravity. The quantities $\alpha$, $g$, $\nu$, $\kappa_T$, $T_{0z}$, and $T_z^{\text{ad}}$ are all assumed to be constant.

We then non-dimensionalize these equations and boundary conditions assuming that the unit length is L, the unit time is
$S_C^{-1}$, and the unit temperature is $L(T_{0z} - T_z^{ad})$, to get

$$ 
\begin{aligned}
\frac{\partial \tilde{u}}{\partial t} + \tilde{u} \cdot \nabla \tilde{u} = & - \nabla \tilde{p} + \text{Ri}_C \tilde{T} e_z + \frac{1}{\text{Re}_C} \nabla^2 \tilde{u}, \\
\nabla \cdot \tilde{u} = & 0, \\
\frac{\partial \tilde{T}}{\partial t} + \tilde{u} \cdot \nabla \tilde{T} + \tilde{w} = & \frac{1}{\text{Pe}_C} \nabla^2 \tilde{T},
\end{aligned}
\quad (15) $$

and

$$ 
\begin{aligned}
\tilde{v} = \tilde{w} = \tilde{T} = & 0 \text{ at } z = \pm 1/2, \\
\tilde{u} = & \pm 1/2 \text{ at } z = \pm 1/2,
\end{aligned}
\quad (16) $$

where pressure, the spatio-temporal coordinates, and the differential operators have also been implicitly non-dimensionalized, and where

$$ 
\begin{aligned}
\text{Re}_C &= \frac{S_C L^2}{\nu}, \quad \text{Pe}_C = \frac{S_C L^2}{\kappa_T}, \\
\text{and } \text{Ri}_C &= \frac{\alpha g (T_{0z} - T_z^{ad})}{S_C^2} = \frac{N^2}{S_C^2},
\end{aligned}
\quad (17) $$

are the Reynolds, Péclet, and bulk Richardson numbers based on the laminar shear $S_C$.

The asymptotic low Péclet number approximation (Lignières 1999) of the equations listed in (15), which we shall use in this paper, is given by

$$ 
\begin{aligned}
\frac{\partial \tilde{u}}{\partial t} + \tilde{u} \cdot \nabla \tilde{u} = & - \nabla \tilde{p} + \text{Ri}_C \tilde{T} e_z + \frac{1}{\text{Re}_C} \nabla^2 \tilde{u}, \\
\nabla \cdot \tilde{u} = & 0, \\
\tilde{w} = & \frac{1}{\text{Pe}_C} \nabla^2 \tilde{T}.
\end{aligned}
\quad (18) $$

These equations are referred to as the LPN equations hereafter. Finally, we also follow the evolution of a non-dimensional passive scalar field $\tilde{c}$ (representing for example the concentration of a given chemical species), according to

$$ \frac{\partial \tilde{c}}{\partial t} + \tilde{u} \cdot \nabla \tilde{c} = \frac{\tau_c}{\text{Pe}_C} \nabla^2 \tilde{c}, \quad (19) $$

where $\tau_c = \kappa_c / \kappa_T$ is the ratio of the microscopic compositional diffusivity to the thermal diffusivity, and where $\tilde{c}$ satisfies the following Dirichlet conditions:

$$ \tilde{c} = \pm 1/2 \text{ at } z = \pm 1/2. \quad (20) $$

In all that follows, we take $\tau_c = 10^{-3}$, which is much larger than realistic values of the diffusivity ratio in stellar interiors (where it is closer to $10^{-9}-10^{-7}$), but the smallest value we can realistically achieve in fully resolved DNSs.

Note that the parameters $\text{Re}_C$, $\text{Pe}_C$ and $\text{Ri}_C$ are different from $\text{Pe}_L$, $\text{Re}_L$ and $J$ defined and used by Prat & Lignières (2013, 2014) and Prat et al. (2016), since the former are based on the laminar shearing rate $S_C$ while the latter are based on the actual shearing rate $S$. As such, one cannot directly compare the results of our simulations at fixed $\text{Re}_C$ and $\text{Ri}_C\text{Pe}_C$ with those of their simulations at the same values of $\text{Re}_L$ and $J\text{Pe}_L$, respectively. However, once the actual value of the non-dimensional shear $\tilde{S} = S/S_C$ in the bulk of the fluid between the plates is measured in our plane Couette simulations, we can calculate what the equivalent values of $\text{Pe}_L$, $\text{Re}_L$ and $J$ would be

using

$$ 
\begin{aligned}
\text{Pe}_L &= \frac{SL^2}{\kappa_T} = \tilde{S}\text{Pe}_C, \quad \text{Re}_L = \frac{SL^2}{\nu} = \tilde{S}\text{Re}_C \\
\text{and } J &= \frac{N^2}{S^2} = \frac{\text{Ri}_C}{\tilde{S}^2},
\end{aligned}
\quad (21) $$

and use these to compare our results to those of Prat and collaborators.

### 3. Numerical Simulations

We solve the set of equations given in (18) in a cubic domain of non-dimensional side-length $L=1$ using a code developed by Verhoeven & Stellmach (2014; see also Verhoeven et al. 2015), modified to include the moving boundaries and to solve the LPN Equations (18) together with (19). The code is doubly periodic and therefore spectral in the horizontal directions, and uses fourth-order finite differences in the vertical direction. We verified that the LPN equations are a good representation of the dynamics of the full Equations (15) in the limit of low $\text{Pe}_C$, as in Garaud & Kulenthirarajah (2016). Table 1 summarizes our results.

The code being originally written for studies of Rayleigh-Bénard convection, the numerical input parameters are the Rayleigh number $\text{Ra} = \text{Ri}_C\text{Pe}_C\text{Re}_C$, Pr = $\text{Pe}_C/\text{Re}_C$ and $\text{Re}_C$. As a result, our parameter sweeps are at constant Ra and constant $\text{Re}_C$ rather than at constant $\text{Ri}_C\text{Pe}_C$ and constant $\text{Re}_C$. Nevertheless, we can still use the data to make useful comparisons between runs. Finally, in all that follows, we use the notations

$$ \bar{q}(z, t) = \frac{1}{L^2} \iint q(x, y, z, t) dxdy, \quad (22) $$

$$ \langle q \rangle (t) = \frac{1}{L^3} \iiint q(x, y, z, t) dxdydz, \quad (23) $$

$$ (q)_t(x, y, z) = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} q(x, y, z, t) dt, \quad (24) $$

to imply horizontal-, volume-, and time-averages, respectively. In the expression for the time average, the times $t_1$ and $t_2$ are selected to bracket the longest possible time-interval available in the simulations during which the shear flow has achieved a statistically stationary state.

#### 3.1. Sample Results

Figure 2 shows volume-rendered snapshots of the vertical velocity field in two simulations at $\text{Re}_C = 1.2 \times 10^5$ and $\text{Ri}_C\text{Pe}_C = 0.0833$ and 83.3, respectively. Note how the typical vertical coherence of the eddies is much larger in the more weakly stratified case (at lower $\text{Ri}_C\text{Pe}_C$) than in the more strongly stratified one, even though it retains significant finescale structure owing to the large Reynolds number of the simulations.

In order to analyze the system more quantitatively, Figure 3 shows the horizontally averaged and time-averaged mean flow $(\tilde{u})_t$ for two sets of simulations: varying $\text{Ri}_C\text{Pe}_C$ at constant $\text{Re}_C$ (left) and varying $\text{Re}_C$ at roughly constant $\text{Ri}_C\text{Pe}_C$ (right). Similarly, Figure 4 shows the horizontally averaged and time-averaged turbulent kinetic energy,

$$ E_{\text{turb}}(z) = \frac{1}{2} (\overline{\tilde{u}^2} - \tilde{u}^2 + \overline{\tilde{v}^2} + \overline{\tilde{w}^2})_t, \quad (25) $$
<table><caption>**Table 1**
Summary of the Main Results for All the Runs</caption><thead><tr><th>Re<sub>C</sub>/10<sup>4</sup></th><th>Ri<sub>C</sub>Pe<sub>C</sub></th><th>N<sub>x</sub>, N<sub>y</sub>, N<sub>z</sub></th><th>S&#773;</th><th>G&#773;</th><th>l<sub>e</sub></th><th>⟨w&#773;c&#773;⟩<sub>t</sub> × 10<sup>4</sup></th><th>⟨w&#773;u&#773;⟩<sub>t</sub> × 10<sup>4</sup></th></tr></thead><tbody><tr><td>4</td><td>2.5 × 10<sup>-4</sup></td><td>96<sup>2</sup> × 128</td><td>0.235 ± 0.003</td><td>0.22 ± 0.002</td><td>0.178</td><td>7.12 ± 1.4</td><td>6.12 ± 2.8</td></tr><tr><td>4</td><td>2.5 × 10<sup>-3</sup></td><td>96<sup>2</sup> × 128</td><td>0.221 ± 0.003</td><td>0.22 ± 0.002</td><td>0.162</td><td>7.12 ± 1.4</td><td>6.12 ± 2.8</td></tr><tr><td>4</td><td>2.5 × 10<sup>-2</sup></td><td>96<sup>2</sup> × 128</td><td>0.246 ± 0.003</td><td>0.23 ± 0.002</td><td>0.178</td><td>7.08 ± 1.4</td><td>6.12 ± 2.8</td></tr><tr><td>4</td><td>2.5 × 10<sup>-1</sup></td><td>96<sup>2</sup> × 128</td><td>0.25 ± 0.003</td><td>0.23 ± 0.002</td><td>0.162</td><td>5.88 ± 1.2</td><td>6.06 ± 2.6</td></tr><tr><td>4</td><td>2.5</td><td>96<sup>2</sup> × 128</td><td>0.345 ± 0.003</td><td>0.29 ± 0.002</td><td>0.116</td><td>5.22 ± 0.8</td><td>5.18 ± 1.7</td></tr><tr><td>4</td><td>8.25</td><td>96<sup>2</sup> × 128</td><td>0.461 ± 0.002</td><td>0.38 ± 0.0015</td><td>0.093</td><td>4.32 ± 0.50</td><td>4.18 ± 1.1</td></tr><tr><td>4</td><td>25</td><td>96<sup>2</sup> × 128</td><td>0.54 ± 0.001</td><td>0.48 ± 0.001</td><td>0.077</td><td>3.37 ± 0.27</td><td>3.25 ± 0.6</td></tr><tr><td>4</td><td>82.5</td><td>96<sup>2</sup> × 128</td><td>0.68 ± 0.001</td><td>0.67 ± 0.001</td><td>0.062</td><td>1.99 ± 0.14</td><td>2.01 ± 0.3</td></tr><tr><td>4</td><td>165</td><td>96<sup>2</sup> × 128</td><td>0.79 ± 0.001</td><td>0.82 ± 0.001</td><td>0.054</td><td>1.03 ± 0.10</td><td>1.10 ± 0.3</td></tr><tr><td>6</td><td>1.66 × 10<sup>-2</sup></td><td>144<sup>2</sup> × 192</td><td>0.224 ± 0.004</td><td>0.21 ± 0.002</td><td>0.27</td><td>5.50 ± 1.0</td><td>5.58 ± 2.3</td></tr><tr><td>6</td><td>1.66 × 10<sup>-1</sup></td><td>144<sup>2</sup> × 192</td><td>0.224 ± 0.006</td><td>0.24 ± 0.004</td><td>0.15</td><td>5.38 ± 0.9</td><td>5.55 ± 2.1</td></tr><tr><td>6</td><td>1.66</td><td>144<sup>2</sup> × 192</td><td>0.313 ± 0.004</td><td>0.28 ± 0.004</td><td>0.10</td><td>4.91 ± 0.8</td><td>4.89 ± 1.2</td></tr><tr><td>6</td><td>16.6</td><td>144<sup>2</sup> × 192</td><td>0.505 ± 0.002</td><td>0.44 ± 0.002</td><td>0.078</td><td>3.41 ± 0.3</td><td>3.2 ± 0.6</td></tr><tr><td>6</td><td>50</td><td>144<sup>2</sup> × 192</td><td>0.602 ± 0.001</td><td>0.55 ± 0.001</td><td>0.067</td><td>2.48 ± 0.1</td><td>2.3 ± 0.3</td></tr><tr><td>6</td><td>166</td><td>144<sup>2</sup> × 192</td><td>0.746 ± 0.001</td><td>0.72 ± 0.0008</td><td>0.047</td><td>1.25 ± 0.07</td><td>1.3 ± 0.2</td></tr><tr><td>9</td><td>1.11 × 10<sup>-3</sup></td><td>192<sup>2</sup> × 256</td><td>0.225 ± 0.004</td><td>0.21 ± 0.003</td><td>0.22</td><td>5.28 ± 1.0</td><td>5.32 ± 1.9</td></tr><tr><td>9</td><td>1.11 × 10<sup>-2</sup></td><td>192<sup>2</sup> × 256</td><td>0.191 ± 0.003</td><td>0.20 ± 0.002</td><td>0.22</td><td>5.39 ± 0.9</td><td>5.36 ± 2.4</td></tr><tr><td>9</td><td>1.11 × 10<sup>-1</sup></td><td>192<sup>2</sup> × 256</td><td>0.214 ± 0.002</td><td>0.21 ± 0.002</td><td>0.15</td><td>5.21 ± 1.0</td><td>5.21 ± 2.1</td></tr><tr><td>9</td><td>1.11</td><td>192<sup>2</sup> × 256</td><td>0.297 ± 0.002</td><td>0.27 ± 0.0015</td><td>0.12</td><td>4.77 ± 0.8</td><td>4.71 ± 1.7</td></tr><tr><td>9</td><td>11.1</td><td>192<sup>2</sup> × 256</td><td>0.480 ± 0.001</td><td>0.40 ± 0.005</td><td>0.085</td><td>3.38 ± 0.3</td><td>3.16 ± 0.6</td></tr><tr><td>9</td><td>111</td><td>192<sup>2</sup> × 256</td><td>0.686 ± 0.0003</td><td>0.64 ± 0.002</td><td>0.047</td><td>1.67 ± 0.07</td><td>1.56 ± 0.2</td></tr><tr><td>9</td><td>333</td><td>192<sup>2</sup> × 256</td><td>0.812 ± 0.0002</td><td>0.80 ± 0.0015</td><td>0.035</td><td>0.84 ± 0.03</td><td>0.74 ± 0.07</td></tr><tr><td>12</td><td>8.33 × 10<sup>-3</sup></td><td>192<sup>2</sup> × 256</td><td>0.204 ± 0.004</td><td>0.20 ± 0.002</td><td>0.19</td><td>4.98 ± 0.9</td><td>4.93 ± 1.9</td></tr><tr><td>12</td><td>8.33 × 10<sup>-2</sup></td><td>192<sup>2</sup> × 256</td><td>0.202 ± 0.003</td><td>0.20 ± 0.002</td><td>0.20</td><td>5.10 ± 1.1</td><td>5.10 ± 2.3</td></tr><tr><td>12</td><td>8.33 × 10<sup>-1</sup></td><td>192<sup>2</sup> × 256</td><td>0.241 ± 0.003</td><td>0.23 ± 0.0015</td><td>0.15</td><td>4.81 ± 0.9</td><td>4.75 ± 1.8</td></tr><tr><td>12</td><td>8.33</td><td>192<sup>2</sup> × 256</td><td>0.464 ± 0.001</td><td>0.38 ± 0.0008</td><td>0.097</td><td>3.50 ± 0.3</td><td>3.21 ± 0.7</td></tr><tr><td>12</td><td>83.3</td><td>192<sup>2</sup> × 256</td><td>0.655 ± 0.0003</td><td>0.59 ± 0.0004</td><td>0.05</td><td>1.83 ± 0.07</td><td>1.65 ± 0.2</td></tr><tr><td>12</td><td>250</td><td>192<sup>2</sup> × 256</td><td>0.762 ± 0.001</td><td>0.74 ± 0.0003</td><td>0.039</td><td>1.05 ± 0.02</td><td>1.00 ± 0.09</td></tr><tr><td>12</td><td>500</td><td>192<sup>2</sup> × 256</td><td>0.850 ± 0.0001</td><td>0.84 ± 0.0002</td><td>0.031</td><td>0.51 ± 0.01</td><td>0.52 ± 0.05</td></tr></tbody></table>

**Note.** The first column reports Re<sub>C</sub>, and the second column reports Ri<sub>C</sub>Pe<sub>C</sub>. The third column is the number of actual or equivalent mesh points in each direction. The fourth column is S&#773;, the amplitude of the mean non-dimensional shear in the middle of the domain (see Section 3.1), and its estimated error, the fifth column is G&#773;, the amplitude of the mean non-dimensional concentration gradient in the middle of the domain (see Section 3.2), and its estimated error, the sixth column is the eddy lengthscale l<sub>e</sub> (see Section 4.1), the seventh column is the time-averaged and volume-averaged turbulent compositional flux and its rms fluctuations, and the eighth column is the time-averaged and volume-averaged Reynolds stress and its rms fluctuations.

as a function of z in the same sets of simulations. Both figures show that the computational domain is divided into a bulk region where the time-averaged shear and turbulent kinetic energy are roughly constant, and two thin boundary layers on either side where the mean flow velocity is forced to match the wall velocity, and where the turbulent kinetic energy first grows somewhat then drops to zero to match the boundary conditions. The tendency for more strongly stratified turbulent systems to have large regions of constant shear is interestingly reminiscent of the numerical results of Garaud & Kulenthirarajah (2016), and will be discussed in more detail in a forthcoming publication.

As the stratification (measured via Ri<sub>C</sub>Pe<sub>C</sub>) increases while holding Re<sub>C</sub> constant, the shear in the bulk of the domain also increases while the turbulent kinetic energy decreases. Meanwhile for fixed Ri<sub>C</sub>Pe<sub>C</sub> ≃ 8.5, increasing the Reynolds number decreases the size of the viscous boundary layer but does not seem to affect the strength of the bulk shear.

Figure 5 shows the amplitude of the non-dimensional shear in the middle of the domain S&#773;, as a function of the input Richardson-Péclet number Ri<sub>C</sub>Pe<sub>C</sub> and of Re<sub>C</sub>. To compute S&#773; we take the time- and horizontally averaged profiles (u&#773;)<sub>t</sub>(z) shown in Figure 3, and fit each of them to the linear function f(z) = -S&#773;z for z ∈ [-0.1, 0.1]. The shear values thus extracted, together with their estimated error, are also reported in Table 1. Note that, with

this definition, S&#773; is equal to minus the gradient of (u&#773;)<sub>t</sub>(z). We see that there are actually two distinct regimes, for low and high stratification, respectively. For very low stratification (as measured by Ri<sub>C</sub>Pe<sub>C</sub> ≪ 1), S&#773; is more-or-less independent of both Re<sub>C</sub> and Ri<sub>C</sub>Pe<sub>C</sub>. As Ri<sub>C</sub>Pe<sub>C</sub> increases beyond one, the effects of stratification clearly become important and a new scaling law emerges with S&#773; ≃ 0.32(Ri<sub>C</sub>Pe<sub>C</sub>)<sup>1/6</sup>. The origin of this law remains to be determined, but likely depends on the structure of the boundary layers near the walls, and is therefore probably model-specific (i.e., in a different model that uses a different type of forcing, the relationship between the midpoint shear and the input Richardson-Péclet number is likely to be different). More importantly, however, we see that it is independent of the Reynolds number, which corroborates our observations of the bulk shear S&#773; in Figure 3. Finally, as Ri<sub>C</sub>Pe<sub>C</sub> continues to increase, S&#773; → 1, which corresponds to the laminar shear solution.

## 3.2. Mixing

One of the goals of this investigation, as discussed in Section 1, is to quantify mixing by shear instabilities, both of a passive tracer and of momentum. A simple model for mixing consists in assuming that the turbulent fluxes are proportional to the local gradient of the mixed quantity, the coefficient of proportionality between the two being called turbulent diffusivity (D<sub>turb</sub>) or turbulent viscosity (ν<sub>turb</sub>). In order to
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Volume-rendered snapshots of the vertical velocity field $\tilde{w}$ in two simulations at $Re_c = 1.2 \times 10^5$. Top: $Ri_c Pe_c = 0.0833$. Bottom: $Ri_c Pe_c = 83.3$.</figcaption></figure>

obtain good statistics, we take a time average and volume average of the turbulent fluxes, and let

$$ \tilde{D}_{\text{turb}} = \frac{\langle \tilde{w} \tilde{c} \rangle_t}{\tilde{G}}, \qquad (26) $$

$$ \tilde{v}_{\text{turb}} = \frac{\langle \tilde{u} \tilde{w} \rangle_t}{\tilde{S}}, \qquad (27) $$

where $\tilde{D}_{\text{turb}} = D_{\text{turb}}/\Delta UL$ is the non-dimensional turbulent diffusivity (and similarly for $\tilde{v}_{\text{turb}}$), and $\tilde{G}$ is minus the gradient of composition in the middle of the domain, measured exactly as $\tilde{S}$ (see Section 3.1) but using the time- and horizontally averaged compositional profile $(\tilde{c})_t(z)$ instead of $(\tilde{u})_t(z)$. Note that the absence of the usual negative sign in these definitions comes from the fact that we have defined $\tilde{G}$ and $\tilde{S}$ to be equal to minus the gradients of $(\tilde{c})_t(z)$ and $(\tilde{u})_t(z)$, respectively.

Figures 6(a) and (b) use the calculated values of $\tilde{D}_{\text{turb}}$ and $\tilde{v}_{\text{turb}}$ to reproduce Figures 3 and 8 of Prat et al. (2016), albeit with the results of our own simulations. Figure 6(a) shows the non-dimensional quantity $D_{\text{turb}}J/\kappa_T = \tilde{D}_{\text{turb}}JPe_c$ as a function of $JPr$, where $J$ is calculated from $\tilde{S}$ as given in (21). We also

show, as a solid line, the parametric model of Equation (12), which was proposed by Prat et al. (2016) to fit their data. We see that, for larger values of $JPr$, our data falls on their parametric curve. In particular, we find as they do that $D_{\text{turb}}J/\kappa_T$ becomes independent of the Reynolds number for large $JPr$ and rapidly drops to zero as $JPr \to (JPr)_c \simeq 0.007$. This confirms the findings of Prat et al. (2016) that diffusive stratified shear instabilities can only occur as long as $JPr < (JPr)_c \simeq 0.007$. Note that this result, together with Equation (5), provides an important experimental constraint on the constants $(JPe)_c$ and $Re_c$.

We also find, as they do, that the data deviates significantly from their parametric model in the limit of low $JPr$, and that the amplitude of the deviations increases with decreasing Reynolds number. In fact, we have ran a few cases with much lower values of $JPr$ than they did (their lowest value of $JPr$ was around $10^{-4}$, while ours is just below $10^{-6}$), and find that this trend continues as $JPr$ decreases. This confirms that their parametric model (and Zahn's model) fares poorly in the limit of low stratification, a result we expected from the discussion of Section 1.2. We revisit the model with the goal of improving it in this limit in Section 4.

Figure 6(b) shows $\nu_{\text{turb}}/\nu = \tilde{v}_{\text{turb}}Re_c$ against $D_{\text{turb}}/\nu = \tilde{D}_{\text{turb}}Re_c$, and illustrates that $\nu_{\text{turb}}$ is indeed proportional to $D_{\text{turb}}$, with a constant of proportionality that varies between 0.8 and 1, as in Prat et al. (2016). Figures 6(a) and (b), and their resemblance with Figures 3 and 8 of Prat et al. (2016), confirm that the properties of the turbulent mixing coefficients $D_{\text{turb}}$ and $\nu_{\text{turb}}$ depend more on the properties of the local shear in the more strongly stratified limit than on the manner in which the shear is forced (i.e., body-forcing versus forcing by surface stresses), the same conclusion reached by Prat et al. (2016) when comparing two different methods of body-forcing (shearing sheet versus explicit forcing). This is reassuring and suggests that, at least in these cases, a local model for turbulent mixing by shear instabilities that only depends on the local gradient Richardson number $J$, together with the local microscopic diffusion coefficients $\nu$ and $\kappa_T$, is appropriate.

## 4. Analysis and Discussion of the Results

We now investigate our data in comparison with Zahn's model (Zahn 1974, 1992) in more depth, and attempt to resolve the outstanding issues described in Section 1.

### 4.1. Comparison of the Various Model Lengthscales

As discussed in Sections 1.1 and 1.2, Zahn's theory relies on the idea that turbulent eddies in diffusive stratified shear flows are just marginally unstable according to the criterion (5). This defines the Zahn scale, $l_Z$, to be (9):

$$ l_Z = \left( \frac{(JPe)_c}{JPe_L} \right)^{1/2} L, \qquad (28) $$

which can then be used in (7) to estimate $D_{\text{turb}}$ (and similarly $\nu_{\text{turb}}$ since they are dimensionally equivalent). As implied by Zahn (1974), however, for turbulence to exist $l_Z$ has to be larger than the scale below which viscosity becomes important, which we call the viscous scale $l_\nu$. The value of $l_\nu$ can be derived from the condition that the eddy-scale Reynolds number is just equal
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left: horizontally averaged and time-averaged mean flow ($\tilde{a}_i$) as a function of $z$, for various simulations at $Re_c = 1.2 \times 10^5$. Right: Same quantity, but this time varying $Re_c$ and holding $Re_i Re_c$ roughly constant. The runs shown have $Re_i Re_c = 8.25, 11.1$ and $8.33$, respectively, as $Re_c$ increases.</figcaption></figure>

to $Re_c$ (see Section 1.1):

$$ \frac{SI_{\nu}^{2}}{\nu} = Re_{c} \Leftrightarrow l_{\nu} = \left(\frac{Re_{c}\nu}{S}\right)^{1/2} = \left(\frac{Re_{c}}{Re_{L}}\right)^{1/2} L. \quad (29) $$

It is easy to verify that the condition $l_Z = l_\nu$ is equivalent to Zahn's stability criterion (5). When this happens, turbulence disappears and $D_{turb}$ must drop to zero instead of being given by (10).

In the limit of weak stratification on the other hand, it is easy to see from (28) that the Zahn scale increases as stratification (or equivalently, $J$) decreases. In a numerical experiment with finite domain size $L$, or in real stellar conditions where the shear layer itself has a finite width $L$, $l_Z$ must eventually approach $L$ as $J$ decreases, at which point the vertical eddy size becomes limited by the domain size instead of being determined by $l_Z$. When this happens, we expect Zahn's mixing prescription to fail. Instead, $D_{turb}$ presumably becomes proportional to $SL^2$ instead of $SI_Z^2$.

To test these ideas, we must compute the height of turbulent eddies $l_e$ and compare it to $L$, $l_Z$, and $l_\nu$. From the discussion above, we predict that (1) $l_e$ should be equal to $L$ (or to a fraction of $L$, as proposed by Prat et al. 2016) in very weakly stratified systems where the eddy size is limited by the system size; (2) $l_e$ should be equal to $l_Z$ in an intermediate regime; and (3) turbulence should disappear in strongly stratified systems when $l_e$ approaches $l_\nu$. The computation of $l_e$ is non-trivial, however. Garaud & Kulenthirarajah (2016) and Prat et al. (2016), who both used models in which perturbations are triply periodic, estimated the eddy scale using a weighted average of the energy spectrum, as in

$$ l_e = \left[ \frac{\sum_{\tilde{k}_s} \sum_{\tilde{k}_v} \sum_{\tilde{k}_z \neq 0} E_{\tilde{w}}(\tilde{k}, t) \tilde{k}_z^{-1}}{\sum_{\tilde{k}_s} \sum_{\tilde{k}_v} \sum_{\tilde{k}_z \neq 0} E_{\tilde{w}}(\tilde{k}, t)} \right]_t, \quad (30) $$

or variants thereof (see Garaud & Kulenthirarajah 2016; Prat et al. 2016, for detail). In this expression, $E_{\tilde{w}}(\tilde{k}, t)$ is the energy in the Fourier component of wavenumber $\tilde{k} = (\tilde{k}_x, \tilde{k}_y, \tilde{k}_z)$ of the vertical velocity field.

In a system with impermeable boundaries such as ours, however, this definition cannot be used as is. In any case, as discussed in Section 1.2, both Garaud & Kulenthirarajah (2016) and Prat et al. (2016) found that using this estimate of the eddy scale leads to inconsistencies in the model, namely that while $D_{turb}$ is approximately predicted by Zahn's model for intermediate values of $J$, it is not proportional to $SI_e^2$ even though this is one of the basic assumptions of Zahn's theory. This statement is equivalent to saying that $D_{turb} = \beta SI_Z^2$ even though $D_{turb} \neq \beta SI_e^2$, implying that $l_e \neq l_Z$. A possible solution of this conundrum is that the turbulent eddies are indeed of size $l_e = l_Z$, but that the definition given in (30) is not a good estimate for $l_e$. This notion is supported by the observation that while the strongly stratified runs are dominated by small-scale eddies (as expected), the more weakly stratified runs, especially at large Reynolds numbers, contain both small and large scales (see Figures 2 and 7). The Zahn scale, presumably, would correspond to the larger one, but the formula given in (30) is strongly weighted toward the smaller scales and would pick these up rather than $l_Z$.

For these reasons, we now propose an alternative way of measuring the eddy scale, using autocorrelation functions of the turbulent flow field. We define

$$ a_{\nu}(l) = \left[ \frac{1}{L-l} \int_{0}^{L-l} \tilde{v}(0, 0, z, t) \tilde{v}(0, 0, z+l, t) dz \right]_t, \quad (31) $$

$$ a_{w}(l) = \left[ \frac{1}{L-l} \int_{0}^{L-l} \tilde{w}(0, 0, z, t) \tilde{w}(0, 0, z+l, t) dz \right]_t, \quad (32) $$

as the mean vertical autocorrelation functions of $\tilde{v}$ and $\tilde{w}$, respectively. The streamwise velocity $\tilde{u}$ is less useful for computing the eddy scale, as it is dominated by the effects of the mean flow even if the latter is subtracted. Note that we integrate over flow profiles taken at $x = 0$ and $y = 0$, rather than over the entire flow field, because those are stored at very regular time intervals in our simulations while the whole flow
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left: Horizontally averaged and time-averaged turbulent kinetic energy $E_{\text{turb}}$ as a function of $z$, for the same simulations as in Figure 3. Right: Same quantity, but this time varying $Re_c$ and holding $Ri_c Pe_c$ roughly constant. The runs shown have $Ri_c Pe_c = 8.25$, 11.1 and 8.33, respectively, as $Re_c$ increases.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Variation of the non-dimensional bulk shearing rate $\tilde{S}$ with $Ri_c Pe_c$. Various symbols correspond to various Reynolds numbers, as seen in the legend. A laminar flow would have $\tilde{S} = 1$. For low $Ri_c Pe_c$, $\tilde{S}$ shows some variation but is overall fairly constant and roughly equal to 0.2. The thick yellow line shows the high $Ri_c Pe_c$ scaling law.</figcaption></figure>

field is only stored once in a while. Since the system is horizontally periodic, and statistically stationary between $t_1$ and $t_2$, we expect that these individual profiles have similar statistics to those of the entire flow; this is actually verified to be true in the Appendix.

There are many possible ways of constructing an eddy lengthscale from these autocorrelation functions. Two of them are explored in the Appendix, and their pros and cons are discussed. We have eventually settled on defining $l_e$ as the first zero of $a_v(l)$:

$$a_v(l_e) = 0. \qquad (33)$$

Figure 8 shows $l_e/L$ thus measured as a function of $JPe_L$ for all the available runs. Generally speaking, we see that $l_e$ decreases with increasing stratification, as expected. In the limit of very large stratification ($JPe_L > 100$), all the data for $l_e$ asymptotes to the same line in this log-log graph. Fitting the data, we find that $l_e/L \simeq (0.5/JPe_L)^{1/2}$, which in turn shows that $l_e$ is proportional to the Zahn scale, since $l_Z/L = ((JPe)_c/JPe_L)^{1/2}$. In fact, if we set the two equal to one another in this limit, in

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** These figures reproduce Figures 3 and 8 of Prat et al. (2016). Top: variation of $D_{\text{turb}}J/\kappa_T = D_{\text{turb}}JPe_c$ against $JPr$, for various Reynolds numbers. Also shown is the parametric model of Prat et al. (2016; see Equation (12)), which fits the data well for the more strongly stratified simulations. Bottom: variation of $\nu_{\text{turb}}/\nu = \dot{v}_{\text{turb}} Re_c$ with $D_{\text{turb}}/\nu = D_{\text{turb}} Re_c$. Note that this is shown on a log-log plot. Also shown are the relationships $\nu_{\text{turb}} = 0.8D_{\text{turb}}$ and $\nu_{\text{turb}} = D_{\text{turb}}$. The former fits the data at intermediate values of the stratification, while the latter fits the data well for the more weakly and more strongly stratified simulations.</figcaption></figure>

line with Zahn's model (Zahn 1992), we can calibrate the unknown constant $(JPe)_c$ to be

$$(JPe)_c \simeq 0.5. \qquad (34)$$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Profiles of $\tilde{w}(0, 0, z)$ at randomly selected times during the quasi-stationary state, for $Re_c = 1.2 \times 10^5$, in the two simulations illustrated in Figure 2. Note how both profiles contain features on small scales, but only the more weakly stratified one contains features on larger-scales too.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between the various theoretical and experimental lengthscales $l_Z$ (solid black line), $l_\nu$ (solid colored lines), and $l_e$ (data points), respectively. The constant ($JPe_c$) in the theoretical Zahn scale (see Equation (28)) has been calibrated to fit the data in the limit of large stratification. That fit then uniquely defines the constant $Re_c$, which is needed to compute $l_\nu$. The colors of the lines and symbols are the same for $l_\nu$ and $l_e$ that share the same value of $Re_c$.</figcaption></figure>

Figure 8 shows the calibrated Zahn scale in addition to the turbulent eddy scale data.

Having calibrated ($JPe_c$), we can also determine $Re_c$ since we had previously found that $(JPr)_c = (JPe)_c / Re_c \approx 0.007$ (see Equation (5)). This implies

$$Re_c \approx 71. \qquad (35)$$

Finally, now that $Re_c$ is known, we can plot $l_\nu$ in Figure 8. For each set of runs (i.e., for each value of the input Reynolds number $Re_c$), we see that turbulent solutions only exist up to the point where $l_Z = l_\nu$ but not beyond. This is entirely consistent with the combined models of Zahn (1974, 1992).

Looking at progressively more weakly stratified systems, Figure 8 shows that $l_Z$ stops being a good estimate for the eddy scale $l_e$ when $JPe_L$ decreases below ~100. Instead, $l_e$ tends to a (roughly) constant fraction of the total domain height (about 20%) when $JPe_L \to 0$. This confirms our suspicion that the local approximation made by Zahn (1992) must fail at low stratification (see also Prat et al. 2016). As discussed earlier, a

better model in this limit would be one in which, for instance, $D_{turb}$ becomes proportional to $SL^2$ rather than $Sl_Z^2$.

In fact, such a model can be obtained using simple dimensional analysis. Indeed, the limit $J \to 0$ is equivalent to the limit $Ri_C \to 0$ (and therefore $Ri_C Pe_C \to 0$) since the two are related via (21). As a result the non-dimensional $\tilde{D}_{turb}$ should asymptote to a value that only depends on the Reynolds number, since it is the only remaining non-dimensional parameter in that limit. If, furthermore, we recall that $Re_c \gg 1$ in stellar interiors, we can anticipate that $\tilde{D}_{turb}$ should tend to a constant $\tilde{D}_0$ that is independent of $Re_c$ as well.<fnref n="4" /> We see in Figure 9(a), which shows $\tilde{D}_{turb}$ as a function of $JPe_C$, that this is indeed the case even for moderate Reynolds numbers, and estimate that $\tilde{D}_0 \approx 0.025$. Dimensionally speaking, this implies $D_{turb} \approx \tilde{D}_0 S_C L^2$, and since $S \approx 0.2 S_C$ in that limit (see Figure 5), we recover $D_{turb} \propto SL^2$ as expected.

In terms of the quantity $D_{turb}J/\kappa_T = \tilde{D}_{turb}JPe_C$, we then have

$$\lim_{J \to 0} \frac{D_{\text{turb}}J}{\kappa_T} = \tilde{D}_0 J Pe_C = \tilde{D}_0 J Pr Re_C, \qquad (36)$$

which is clearly quite different from the behavior of Zahn’s model or of the parametric fit of Prat et al. (2016) in the same limit (where $D_{turb}J/\kappa_T$ tends to a constant and has no dependence on $Re_c$). In Figure 9(b), we see that our data indeed support the notion that $D_{turb}J/\kappa_T \propto JPrRe_c$ in the limit of small $JPr$.

## 4.2. A New Model

We now set out to propose a new model that correctly captures all three limits: the limit of low $J$, where $D_{turb}J/\kappa_T$ should be a function $JPe_C$ by dimensional analysis, the limit of large $J$, where $D_{turb}J/\kappa_T$ should drop to zero when $JPr \to (JPr)_c$, and the intermediate regime where Zahn’s local model (Zahn 1992) directly applies and $D_{turb}J/\kappa_T$ should tend to a constant. This lead us to propose the following functional form:

$$\tilde{D}_{\text{turb}} = f_1(JPe_c) f_2(JPr), \qquad (37)$$

where the functions $f_1$ and $f_2$ satisfy the following properties: (1) $f_2$ tends to 1 for small $JPr$, (2) $f_2$ tends to 0 as $JPr$ tends to $(JPr)_c \approx 0.007$, (3) $f_1$ tends to about $\tilde{D}_0 \approx 0.025$ as $JPe_C$ tends to 0, and (4) $f_1$ tends to $C(JPe_C)^{-1}$ for large $JPe_C$ (but small $JPr$). The first condition is simply a normalization condition. The second is to satisfy the empirical stability criterion. The third is to match the data in the limit of small $JPe_C$ and the fourth, finally, is to recover Zahn’s model in the intermediate limit where $JPe_C$ is large but $JPr$ is small. After some experimentation, we have found that the following model matches the data quite well, and is reasonably simple:

$$\tilde{D}_{\text{turb}} = \frac{C}{a + JPe_C} \left(1 - \frac{JPr}{(JPr)_c}\right)^b, \quad \text{if } JPr < (JPr)_c, \\ \tilde{D}_{\text{turb}} = 0 \text{ if } JPr > (JPr)_c, \qquad (38)$$

<footnote n="4">⁴ Note that $\tilde{D}_0$ is model-specific, in the sense that using a different model setup would yield a different constant.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of the data and the model for the turbulent diffusivity $D_{\text{turb}}$. The top panel shows this comparison in terms of the non-dimensional diffusivity $\tilde{D}_{\text{turb}}$ vs. $J\text{Pe}_C$, which emphasizes the behavior of the system in the weakly stratified limit. The bottom panel shows the same data and models in terms of $D_{\text{turb}}J/\kappa_T$ vs. $J\text{Pr}$, which emphasizes the behavior of the system in the strongly stratified limit. In both panels the symbols show the data, and the lines of the corresponding color show the model predictions for the same $\text{Re}_C$. Also shown are model predictions for $\text{Re}_C = 10^6$ and $\text{Re}_C = 10^7$. Note that the error bars have been omitted to avoid crowding the plot. The thick black line represents Zahn's basic model (10) with $C = 0.08$, see Section 4.2.</figcaption></figure>

or equivalently

$$ \begin{gathered} \frac{D_{\text{turb}}J}{\kappa_T} = \frac{C}{1 + a(J\text{Pe}_C)^{-1}} \left(1 - \frac{J\text{Pr}}{(J\text{Pr})_c}\right)^b \quad \text{if } J\text{Pr} < (J\text{Pr})_c, \\ \frac{D_{\text{turb}}J}{\kappa_T} = 0 \quad \text{if } J\text{Pr} > (J\text{Pr})_c. \end{gathered} \tag{39} $$

Note that $C/a = \tilde{D}_0 \simeq 0.025$, so the constants $C$ and $a$ are not independent of one another. Fitting the data, we find that using

$$ C \simeq 0.08, \quad a \simeq 30, \quad b \simeq 0.25, \quad \text{and} \quad (J\text{Pr})_c \simeq 0.007 \tag{40} $$

provides a reasonably good overall fit to all of our simulations.<fnref n="5" /> Note that having found $C$, we can then obtain the constant $\beta$ of Zahn's model (see Equation (7)): since $C = \beta(J\text{Pe})_c$ (see Equation (10)) we have

$$ \beta \simeq 0.16. \tag{41} $$

The model is shown in Figure 9, for various values of the Reynolds number corresponding to the numerical simulations available, as well as for much larger Reynolds numbers for

which simulations are not currently possible, namely $\text{Re}_C = 10^6$ and $\text{Re}_C = 10^7$.

In the limit of very high Reynolds numbers ($\text{Re}_C = 10^7$), the model prediction clearly exhibits the three significant regimes discussed above: the very high $J$ cutoff, in which the model only depends on $J\text{Pr}$ and where $D_{\text{turb}}$ drops to zero, the very low $J$ limit where $\tilde{D}_{\text{turb}} \to \tilde{D}_0$ or equivalently $D_{\text{turb}}J/\kappa_T \to \tilde{D}_0J\text{Pr}\text{Re}_C$, and an intermediate region where $J\text{Pe}_C$ is large while $J\text{Pr} = J\text{Pe}_C/\text{Re}_C$ is small, where $D_{\text{turb}}J/\kappa_T \simeq C$, or equivalently $\tilde{D}_{\text{turb}} \simeq C(J\text{Pe}_C)^{-1}$. As the Reynolds number decreases, however, the intermediate region shrinks and the dynamics rapidly go from being domain-size-dominated to viscously dominated. In the $\text{Re}_C = 10^6$ case, for instance, one only barely distinguishes the regime where $D_{\text{turb}}J/\kappa_T \simeq C$. This explains why in our numerical simulations (and similarly in those of Prat et al. 2016), where the largest value of $\text{Re}_C$ is of the order of $10^5$, the intermediate region is effectively absent. Despite this, we see that the new model fits the data very well for all available simulations, both in the strongly stratified limit of course, where it recovers the essence of Zahn's 1974 and 1992 models, but also in the weakly stratified limit where it accounts for the limitation of the eddy scale either by the domain size or by the shear lengthscale.

### 4.3. Discussion of the Proposed Model

The proposed model and its fitted constants, namely those used in (39) as well as the others introduced in Section 1 ($\text{Re}_c$, $J\text{Pe}_c$, and $\beta$), require some discussion.

First, note that $C$, $b$, and $(J\text{Pr})_c$ are tied to the properties of the local model only—the first being from Zahn's original model (10), and the other two being related to the proposed correction in the limit of $J\text{Pr} \to (J\text{Pr})_c$. In this sense, they should be fairly universal and are not expected to depend, say, on the large-scale properties of the shear or on the boundary conditions applied in the model. The constraints on $(J\text{Pr})_c$ are fairly strong given that values close to 0.007 were independently measured in various model setups by Prat et al. (2016), Garaud & Kulenthirarajah (2016), and in this paper. Strict error bars on $(J\text{Pr})_c$ still remain to be estimated from a larger suite of simulations, however. The constraints on $b$ are arguably weaker. In fact, whether the variation of $D_{\text{turb}}J/\kappa_T$ in the limit of $J\text{Pr} \to (J\text{Pr})_c$ is best represented by a power law, as we propose here, or by another function that tends to zero at $J\text{Pr} = (J\text{Pr})_c$, remains to be determined. In the meantime, the proposed power law with $b \simeq 0.25$ seems to be adequate. Finally, it is important to note that the model constraints on $C$ are also quite weak at this point. The constant $C$ uniquely controls the value of $D_{\text{turb}}$ in the intermediate regime where the Zahn scale is much smaller than the domain size, and yet much larger than the viscous scale. Unfortunately, this regime is never achieved with the values of the Reynolds number presently accessible to moderate-scale numerical simulations. Instead, we are forced to fit $C$ under less-than-ideal conditions where it partially contributes to the variation of $D_{\text{turb}}$ in the limits where $J\text{Pr} \to (J\text{Pr})_c$ or $J\text{Pr} \to 0$, together with all the other model constants. Hence uncertainties in estimating $a$, $(J\text{Pr})_c$ and $b$ all affect our estimate for $C$. To address this problem, we plan to run a few simulations at much larger Reynolds number ($\text{Re}_C = O(10^6)$) in the future to better constraint $C$ independently of the other constants.

<footnote n="5">⁵ Note that until more data are obtained in the large Reynolds number limit to improve our constraints on $C$, it is somewhat pointless to attempt to fit these constants within high degrees of accuracy.</footnote>
Second, it is crucial to understand that, by contrast with C, b, and (JPr)$_c$, the constant a is not universal but instead depends on the global properties of the model (such as the shape and amplitude of the mean shear, its lengthscale compared with the domain size, and the manner in which the model is non-dimensionalized). To see this more clearly, consider a thought experiment where we simply non-dimensionalize our plane Couette flow model in two different ways, using, say, the original lengthscale L in one case and L/2 in the other. The prediction for $D_{\text{turb}}J/\kappa_T$ must remain unchanged, since the (dimensional) model setup is the same in both cases. However, since the Péclet number definitions are different (one being Pe$_1$ = Pe$_C$ and the other being Pe$_2$ = Pe$_C$/4), the value of a would have to change by a factor of 4 to compensate (see Equation (39)). The constant a would similarly have to change should one decide to use a different unit time than $S_C^{-1}$. In other words, a is inherently model-dependent.

Hence, if one wishes to apply Equation (39) to predict how much mixing will take place in a system different from the one introduced in Section 2, C, b and (JPr)$_c$ would remain unchanged but a needs to be re-calibrated. To verify that this is indeed the case, and by way of providing an example, we revisit the data from Garaud & Kulenthirarajah (2016) for the sinusoidally body-forced shear flow (see Section 1). We extract the mean shear S at z = π (which is in the middle of the domain and in the middle of the turbulent shear layer in their simulations), as well as the turbulent viscosity<fnref n="6" /> at that position using the technique described in this paper; see Section 3.2. We then compare $\nu_{\text{turb}}J/\kappa_T$ to the model predictions, this time using the formula

$$ \frac{\nu_{\text{turb}}J}{\kappa_T} = \frac{C}{1 + a(J\text{Pe})^{-1}} \left(1 - \frac{J\text{Pr}}{(J\text{Pr})_c}\right)^b \quad \text{for } J\text{Pr} < (J\text{Pr})_c, \qquad (42) $$

where Pe is the input Péclet number based on the laminar flow (see the definition in Garaud & Kulenthirarajah 2016), where $C = 0.08$, $b = 0.25$, $(JPr)_c = 0.007$, but where a has been re-fitted and found to be about 8. Note that we have assumed for simplicity that $\nu_{\text{turb}} \simeq D_{\text{turb}}$. The results are shown in Figure 10. We see that, with this re-fitting, the model adequately captures nearly all the features of the data from the weakly to the strongly stratified limit, giving us confidence that it is indeed the correct approach and that our model is fairly universal (aside from a).

Of course, it is a little disappointing—but unavoidable given its meaning—that there should be a model-dependent constant left. The value of a can be determined from numerical experiments under a specific model setup, but it remains unclear what to select for it in the context of a stellar evolution model where the general shape of the shear, and the forcing mechanisms, vary from star to star. In practice, however, the uncertainty in a should not affect the model predictions too much in the limit of large Reynolds numbers appropriate for stellar interiors. Indeed, as shown in Figure 9, the region of parameter space where the value of a influences $D_{\text{turb}}J/\kappa_T$ is limited to that of weak stratification where $D_{\text{turb}}J/\kappa_T \simeq (C/a)J\text{Pe}_C$, and the size of that region shrinks with increasing Reynolds number. One

<figure><img src="image_12.png" /><figcaption>Figure 10. Comparison of our new model (see Equation (42)) with the data from body-forced sinusoidal shear profiles of Garaud & Kulenthirarajah (2016). The constants $C = 0.08$, $b = 0.25$, $(JPr)_c = 0.007$ have not been changed, but a has been re-fitted and found to be about 8 in this case.</figcaption></figure>

may therefore hope that the values of JPr actually relevant for stars would always fall in either the intermediate regime, or in the strongly stratified regime. In the meantime, for the purpose of picking a value of a to be used in stellar evolution codes, one might as well use something completely generic such as $a = 1$.

Finally, a related item of discussion concerns the model constants $(J\text{Pe})_c$, Re$_c$, and β. By contrast with C, b, and (JPr)$_c$, they are not universal because they depend on the definitions we have used to evaluate the eddy scale $l_e$. To see this, suppose we had instead measured the vertical eddy scale using the autocorrelation function of the vertical velocity field $a_v(l)$ instead of $a_v(l)$ (see Equation (32) and the Appendix for details). Applying the same steps as in Section 4.1, we would have found that $(J\text{Pe})_c$ is roughly equal to 2 instead of 0.5, that Re$_c \simeq 284$, and $\beta = 0.04$. One should therefore always bear the relative arbitrariness of $(J\text{Pe})_c$, Re$_c$, and β in mind. Ultimately, however, this does not pose any problem since our model has been written in such a way as to contain only the constants C, b, and (JPr)$_c$ that are indeed universal, together with a, which was discussed earlier.

## 5. Summary and Future Prospects

In this paper we have analyzed a suite of numerical experiments of diffusive stratified shear instabilities using a plane Couette flow setup, and under both the Boussinesq approximation (Spiegel & Veronis 1960) and the low Péclet number approximation (Lignières 1999). We were able to span a range of Reynolds numbers from $4 \times 10^4$ up to $1.2 \times 10^5$, and a range of Richardson-Péclet number from about $10^{-4}$ to about $10^4$. Our results are very comparable with those of Prat and collaborators (Prat & Lignières 2013, 2014; Prat et al. 2016) even though they used a different model setup, and the favorable comparison can be viewed as a successful validation of the codes and of our selected approaches.

We have analyzed our results in the light of Zahn’s models (Zahn 1974, 1992) and found that (1) the stability criterion of Zahn (1974), given in (5), appropriately describes the upper limit for turbulent mixing in simulations of stratified shear flows, with $(JPr)_c \simeq 0.007$ as initially found by Prat et al. (2016); and (2) the turbulent mixing model of Zahn (1992) given in (10), on the other hand, needs to be augmented to account for the stability cutoff in the strongly stratified limit, and for non-local effects in the weakly stratified limit. For

<footnote n="6">⁶ Note that Garaud & Kulenthirarajah (2016) did not add a passive scalar field for their low Péclet number runs, which is why here we use the turbulent viscosity instead.</footnote>
intermediate values of the stratification, on the other hand, we have provided some evidence that the model is likely appropriate, especially at large Reynolds numbers.

In order to better account for all three regimes (weak, intermediate, and strong stratification), we then proposed a new model for diffusive stratified shear instabilities that extends and unifies the models of Zahn (1974, 1992) and fits the results of our numerical experiments quite well. In this model, the turbulent diffusivity (of a passive tracer) is given by

$$D_{\text{turb}} = \frac{C}{1 + a(J\text{Pe})^{-1}} \left(1 - \frac{J\text{Pr}}{(J\text{Pr})_c}\right)^b \frac{\kappa_T}{J} \quad \text{if } J\text{Pr} < (J\text{Pr})_c, \qquad (43)$$

where Pe = SL²/κT and L is the shear scaleheight. The model constants, which have been fitted to our data, include

$$C \approx 0.08, b \approx 0.25, \text{ and } (J\text{Pr})_c \approx 0.007 \qquad (44)$$

which are universal constants, and a which depends on the model setup considered. Lacking further information, we suggest the use of a = 1 in stellar evolution calculations just for simplicity. While a plays an important role in the model predictions in numerical simulations at moderate Reynolds numbers, it should not really affect the predictions for D<sub>turb</sub>J/κ<sub>T</sub> in the limit of large Reynolds numbers and large J appropriate for stellar interiors (see Section 4.3 for a discussion of this issue). Finally, we have also shown that

$$\nu_{\text{turb}} \approx \gamma D_{\text{turb}}, \qquad (45)$$

where the proportionality constant γ appears to vary between about 0.8 and 1 depending on the regime considered. This result is similar to the findings of Prat et al. (2016). Since this variation could be due to the influence of boundaries in our system, we prefer not to attempt to constrain it further here. Nevertheless, if a simple order-of-magnitude estimate is required, selecting γ = 1 is appropriate.

As discussed in Section 4.3, while (JPr)<sub>c</sub> is reasonably well constrained, significant uncertainties remain on b and C. Uncertainties on b can be reduced by running and analyzing a suite of numerical simulations that systematically explore the region of parameter space where JPr → (JPr)<sub>c</sub>. Uncertainties on C can be reduced by running and analyzing a few numerical simulations at significantly higher values of the Reynolds number than what we have presented here. In both cases, the tasks require significant computational resources, but this will be pursued in the future.

In addition to the uncertainties on the model constants, three additional issues remain to be studied. The first is whether Zahn's models (Zahn 1974, 1992) and their extensions discussed in this paper, which were principally derived and tested under conditions where the shear is constant and where the entire computational domain is turbulent, would also apply in stellar situations where turbulent regions would likely coexist with stable regions, and where the shear may be quite far from being constant. Garaud & Kulenthirarajah (2016) used a sinusoidal body force that generated a non-constant mean shear S(z), but focused on analyzing the global properties of their simulations, namely the volume-averaged turbulent diffusion coefficient, and whether the system as a whole becomes turbulent or not. It would be interesting to revisit their simulations to see whether the new model given in (43) applies

locally as well (i.e., at each vertical position z). The preliminary investigation presented in Figure 10 suggests that it might, at least deep in the middle of the turbulent layer. Whether this continues to be true in all circumstances, especially in turbulent regions bordering stable ones, remains to be determined and will be the subject of Paper 2 in this series.

The second issue concerns the effect of rotation. Zahn's models (Zahn 1974, 1992) are principally invoked in the context of rotational shear, and in fact the expression for the turbulent diffusion coefficient is often written as

$$D_{\text{turb}} = C \frac{\kappa_T}{N^2} \left( r \frac{d\Omega}{dr} \right)^2, \qquad (46)$$

where the shear has been implicitly derived from the rotational shear rdΩ/dr (where Ω(r) is the rotation rate of a shell at radius r in the star). However, rotation was not included in any of the simulations performed to date. This is problematic because rotation can have both stabilizing and destabilizing effects on stratified fluids, depending on the sign of the angular momentum gradient, which could dramatically change the results obtained so far. Hence, future work must include rotation to provide a more comprehensive test of the theories proposed by Zahn, and in this paper.

Finally, the third issue concerns the limit of validity of the model beyond the use of the low Péclet number approximation. In this paper, we have restricted our attention to numerical simulations that use the LPN equations (Lignières 1999). As discussed by Lignières (1999), and confirmed in the numerical simulations of Garaud & Kulenthirarajah (2016), this approximation is valid as long as the turbulent Péclet number of the fluid is smaller than one, but not when it exceeds one. Whether Zahn's models would continue to apply when the LPN approximation is invalid remains to be determined. This is a rather formidable computational problem since it requires having a large Péclet number together with a small Prandtl number, which implies a very large Reynolds number. Nevertheless, it is a crucial problem since in all likelihood a significant fraction of stellar shear layers may actually belong to this region of parameter space rather than to the low Péclet number limit. The dynamics of such high Péclet number shear flows, however, should bear some resemblance with those of stratified shear flows in the geophysical context, so we may be able to use some of the existing literature (Phillips 1969; Peltier & Caulfield 2003; Mahrt 2014) to gain a better understanding of the various processes involved.

P.G. and D.G. gratefully acknowledge funding by NSF AST-1517927. Support was provided for J.V. by the National Aeronautics and Space Administration under grants OPR NNX13AK94G and PGG NNX14AN70G. The simulations were run on the Hyades supercomputer at UCSC, purchased using an NSF MRI grant.

## Appendix

In this appendix we study the autocorrelation functions a<sub>v</sub>(l) and a<sub>w</sub>(l) defined in Equations (31) and (32) and the manner in which they can be used to define a vertical eddy lengthscale.
<figure><img src="image_13.png" /><figcaption>Figure 11. Comparison of $A_v(l, t)$ and $A_w(l, t)$ at 10 individual timesteps (labeled “Whole domain,” shown as thin black lines), with their average (labeled “Whole domain, mean,” shown as the red dashed line), and with $a_v(l)$ and $a_w(l)$ (labeled “Single profiles, mean,” shown in the blue dotted line), for $Re_c = 1.2 \times 10^5$, in the two simulations illustrated in Figure 2 with $Ri_cPe_c = 0.083$ and 83.3, respectively.</figcaption></figure>

We begin by comparing the autocorrelation functions $a_v(l)$ and $a_w(l)$ with

$$A_v(l, t) = \frac{1}{L^2} \iint \left[ \frac{1}{L-l} \int_0^{L-l} \tilde{v}(x, y, z, t) \times \tilde{v}(x, y, z+l, t) dz \right] dxdy, \quad (47)$$

$$A_w(l, t) = \frac{1}{L^2} \iint \left[ \frac{1}{L-l} \int_0^{L-l} \tilde{w}(x, y, z, t) \times \tilde{w}(x, y, z+l, t) dz \right] dxdy, \quad (48)$$

namely the vertical autocorrelation functions at individual times $t$, but integrated over the entire domain.

In theory, the best way to compute the mean vertical autocorrelation of the flow field would be to take the time average of the functions $A_v(l, t)$ and $A_w(l, t)$ once the turbulent shear flow is in a statistically stationary state. In practice, however, the functions $A_v(l, t)$ and $A_w(l, t)$ take significantly longer to compute than $a_v(l)$ and $a_w(l)$, and require storage of the entire flow field at each time $t$, which we only have for a handful of simulations. Computing $a_v(l)$ and $a_w(l)$, on the other hand, only requires storing single vertical profiles at $x = y = 0$ at each time $t$, which is significantly less expensive and was done by default in all the runs. However, one may wonder whether the functions $a_v(l)$ and $a_w(l)$ computed this way are statistically equivalent to the time-averages of $A_v$ and $A_w$ or not.

To test this, we have computed $A_v(l, t)$ and $A_w(l, t)$ at 10 roughly equidistant individual timesteps at which the whole flow field was saved, in the two simulations presented in Figure 2, as well as the algebraic average of $A_v(l, t)$ and $A_w(l, t)$ over these 10 steps. This is of course not the same as taking the time average of $A_v(l, t)$ and $A_w(l, t)$, but is at least an attempt at approximating that time average. We then compare in Figure 11 both individual profiles of $A_v(l, t)$ and $A_w(l, t)$ and their averages to $a_v(l)$ and $a_w(l)$ computed for the same simulations.

Figure 11 reveals a number of things. First, we see that $a_v(l)$ and $a_w(l)$ (dotted blue lines) are indeed a good approximation to the true vertical autocorrelation function of the flow field (approximated here by the dashed red curves) for low and high stratification, at least for values of $l$ smaller than about 0.5. This confirms that using single profiles at $x = y = 0$ to compute the vertical autocorrelation functions of the entire flow field is indeed satisfactory.

Second, looking at $a_v(l)$ and $a_w(l)$ in more detail, we see that they are much more variable, and drop much more slowly with $l$ in the weakly stratified case than in the strongly stratified case. These results are expected from the notion that vertical motion becomes progressively more restricted as stratification increases. Somewhat more surprising is the fact that $a_v(l)$ has a systematically longer “tail” at large $l$ than $a_v(l)$ at the same value of the stratification, and does not necessarily change sign
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Estimate of the vertical eddy size using Equation (49), for all of the simulations presented in Table 1. Note how the data are similar in shape but overall much more scattered than in Figure 8, which uses the definition of (50) instead. Also shown are the re-calibrated Zahn scale (this time with ($JPe_e$) = 2) and corresponding viscous scale $l_v$.</figcaption></figure>

while $a_v(l)$ always does. We believe that this is because we are comparing somewhat different ideas, namely the vertical correlation of vertical flows (in $a_w(l)$) to the vertical correlation of horizontal flows (in $a_v(l)$)—the former is indeed likely to extend further than the latter.

Naively speaking, the most obvious way of defining the eddy lengthscale $l_e$ would be to take the first zero of $a_w(l)$, interpreting this as the scale over which vertical fluid motions change from being upward to downward, or vice-versa. In practice, however, Figure 11 shows that this definition would be problematic since $a_w(l)$ does not always reliably have a zero. One possible way to solve the problem is to define $l_e$ instead as the value of $l$ where $a_w(l)$ first drops below a certain fraction, say 5%, of $a_w(0)$:

$$a_w(l_e) = 0.05a_w(0). \qquad (49)$$

The value of 5% selected here is somewhat arbitrary, in as much as we have found that any value between 5% and 10% yields predictions for $l_e$ that are qualitatively similar. Selecting smaller values picks up too much of the statistical noise, while selecting larger values may not appropriately define an eddy size. As another possibility, having noticed that $a_v(l)$ consistently changes sign at a given lengthscale $l$, we can also define $l_e$ as the

lengthscale for which

$$a_v(l_e) = 0. \qquad (50)$$

The values of $l_e$ measured using (49) are shown in Figure 12, while those using (50) are shown in the main text in Figure 8. In both cases the eddy scale typically decreases with increasing $JPe_L$, and is reasonably independent of the Reynolds number. We can clearly see, however, that the data shown in Figure 12 have much more scatter than in Figure 8, a result that can be attributed to the relatively poor conditioning of Equation (49). Indeed, since $a_w(l)$ has a relatively flat tail, $l_e$ is very sensitive to small variations in $a_w(l)$. For this reason, we ultimately chose to use the definition of $l_e$ based on $a_v(l)$ in this work (see Equation (50)), which does not suffer from this poor conditioning.

## References

Bischoff, T. 2013, in Proc. 2013 Program in Geophysical Fluid Dynamics, ed. C. Cenedese, E. Chassignet, & S. Llewellyn-Smith

Brown, T. M., Christensen-Dalsgaard, J., Dziembowski, W. A., et al. 1989, ApJ, 343, 526

Christensen-Dalsgaard, J., & Schou, J. 1988, in ESA Special Publ. 286, Seismology of the Sun and Sun-Like Stars, ed. E. J. Rolfe (Noordwijk: ESA), 149

Deheuvels, S., et al. 2012, ApJ, 756, 19

Deheuvels, S., et al. 2014, A&A, 564, A27

Dudis, J. J. 1974, JFM, 64, 65

Garaud, P., Gallet, B., & Bischoff, T. 2015, PhFl, 27, 084104

Garaud, P., & Kulenthirarajah, L. 2016, ApJ, 821, 49

Howard, L. N. 1961, JFM, 10, 509

Jones, C. A. 1977, GApFD, 8, 165

Lignières, F. 1999, A&A, 348, 933

Lignières, F., Califano, F., & Mangeney, A. 1999, A&A, 349, 1027

Mahrt, L. 2014, AnRFM, 46, 23

Miles, J. W. 1961, JFM, 10, 496

Peltier, W. R., & Caulfield, C. 2003, AnRFM, 35, 135

Phillips, O. M. 1969, AnRFM, 1, 245

Prat, V., Guilet, J., Viallet, M., & Müller, E. 2016, A&A, 592, A59

Prat, V., & Lignières, F. 2013, A&A, 551, L3

Prat, V., & Lignières, F. 2014, A&A, 566, A110

Richardson, L. F. 1920, RSPSA, 97, 354

Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442

Townsend, A. A. 1958, JFM, 4, 361

Verhoeven, J., & Stellmach, S. 2014, Icar, 237, 143

Verhoeven, J., Wiesehöfer, T., & Stellmach, S. 2015, ApJ, 805, 62

Zahn, J.-P. 1974, in IAU Symp. 59, Stellar Instability and Evolution, ed. P. Ledoux, A. Noels, & A. W. Rodgers (Dordrecht: D. Reidel), 185

Zahn, J.-P. 1992, A&A, 265, 115