# Turbulent Transport by Diffusive Stratified Shear Flows: From Local to Global Models. III. A Closure Model

Logithan Kulenthirarajah, Pascale Garaud

► To cite this version:

Logithan Kulenthirarajah, Pascale Garaud. Turbulent Transport by Diffusive Stratified Shear Flows: From Local to Global Models. III. A Closure Model. The Astrophysical Journal, 2018, 864, 10.3847/1538-4357/aad5e8. insu-03678183

HAL Id: insu-03678183

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

Submitted on 20 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
<figure><img src="image_2.png" /></figure>

# Turbulent Transport by Diffusive Stratified Shear Flows: From Local to Global Models.
## III. A Closure Model

Logithan Kulenthirarajah¹,² and Pascale Garaud³

¹Institut de Recherche en Astrophysique et Planétologie (IRAP), Toulouse, France

²Université Fédérale Toulouse Midi-Pyrénées (UFTMiP), France

³Department of Applied Mathematics, Baskin School of Engineering, University of California Santa Cruz, 1156 High Street, Santa Cruz, CA 95060, USA
Received 2018 March 30; revised 2018 May 15; accepted 2018 May 23; published 2018 September 5

### Abstract

Being able to account for the missing mixing in stellar radiative zones is a key step toward a better understanding of stellar evolution. Zahn argued that thermally diffusive shear-induced turbulence might be responsible for some of this mixing. In Part I and Part II of this series of papers, we showed that Zahn’s mixing model applies when the properties of the turbulence are local. But we also discovered limitations of the model when this locality condition fails, in particular near the edge of a turbulent region. In this paper, we propose a second-order closure model for the transport of momentum and chemical species by shear-induced turbulence in strongly stratified, thermally diffusive environments (the so-called low Péclet number (LPN) limit), which builds upon the work of Garaud & Ogilvie. Comparison against direct numerical simulations (DNSs) shows that the model is able to predict the vertical profiles of the mean flow and of the stress tensor (including the momentum transport) in diffusive shear flows, often with a reasonably good precision, and at least within a factor of order unity in the worst case scenario. The model is sufficiently simple to be implemented in stellar evolution codes, and all the model constants have been calibrated against DNSs. While significant limitations to its use remain (e.g., it can only be used in the slowly rotating LPN limit), we argue that it is more reliable than most of the astrophysical prescriptions that are used in stellar evolution models today.

**Key words:** hydrodynamics – instabilities – stars: evolution – stars: interiors – turbulence

### 1. Introduction

First discussed in the context of stellar astrophysics by Zahn (1974; see also Endal & Sofia 1978; Zahn 1992; Maeder 1995; Maeder & Meynet 2000), diffusive shear instabilities—also called secular shear instabilities—are now considered by most stellar evolution codes (Pinsonneault et al. 1989; Heger et al. 2000; Eggenberger et al. 2008; Paxton et al. 2013) and can be a significant source of mixing in some phases of stellar evolution. This is in contrast with standard shear instabilities (in which perturbations are assumed to be adiabatic), which usually cannot develop in stellar radiative zones owing to the overwhelmingly stabilizing effect of the thermal stratification. Indeed, the well-known energy-based criterion for instability in the adiabatic case is (Richardson 1920)

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

where N is the local Brunt–Väisälä frequency, S is the local shearing rate, and $J_c$ is a constant of order unity. It is rarely satisfied in radiative zones since J is usually much larger than one. Zahn (1974) argued that thermal diffusion can, however, mitigate the stabilizing impact of thermal stratification, to the extent that the relevant criterion for instability should no longer be given by Equation (1) but instead by

$$JPr < (JPr)_c, \qquad (2)$$

where $Pr = \nu/\kappa_T$ is the Prandtl number ($\nu$ is the kinematic viscosity and $\kappa_T$ is the thermal diffusivity) and where $(JPr)_c$ is a number of order $10^{-3}$. Since $Pr \sim 10^{-9}-10^{-6}$ in stellar interiors, Zahn’s criterion would suggest that diffusive shear

instabilities are possible for J up to $10^3-10^6$, which is more in line with the range of observed values in stars.

Zahn (1992) later put forward a simple model for turbulent transport by such instabilities, arguing that the typical scale l of the turbulent eddies is the largest possible scale for which $JPr$ is of order unity, where

$$Pe_l = \frac{Sl^2}{\kappa_T} \qquad (3)$$

is the eddy-scale Péclet number (i.e., the ratio of the advective timescale to the thermal diffusion timescale). That scale will hereafter be referred to as the Zahn scale and is given by

$$J \frac{Sl_z^2}{\kappa_T} = (JPe)_c \Rightarrow l_z = \sqrt{\frac{(JPe)_c \kappa_T}{JS}}, \qquad (4)$$

where $(JPe)_c$ is a constant of order one. An estimate for the turbulent diffusion coefficient can then readily be obtained by simple dimensional analysis (Zahn 1992):

$$D_{\text{turb}} = \beta Sl_z^2 = \beta (JPe)_c \frac{\kappa_T}{J} \equiv C \frac{\kappa_T}{J}, \qquad (5)$$

where $\beta$ and $C$ are thus both constants as well. Applying Zahn’s model for diffusive shear instabilities therefore involves using (5) whenever (2) is satisfied. It is important to note, however, that this model is purely local and is not expected to apply when the Zahn scale becomes commensurate with the shear length scale, or with any other macroscopic system scale, such as the distance to the surface or to the nearest convective zone for instance.

In the past few years, tremendous progress has been made using both theory and Direct Numerical Simulations (DNSs) to
test Zahn's model and to gain a more detailed understanding of the properties of diffusive shear instabilities. As DNSs are limited to modeling a very small region of a star, the Boussinesq approximation for gases (Spiegel & Veronis 1960) is most commonly used, in which one assumes the fluid is in hydrostatic and thermal equilibrium, and where perturbations evolve according to

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = - \frac{\nabla p}{\rho_0} + \alpha g T e_z + \nu \nabla^2 \mathbf{u}, \quad (6) $$

$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T + u_z (T_{0z} - T_{ad,z}) = \kappa_T \nabla^2 T, \quad (7) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (8) $$

where $\mathbf{u} = (u_x, u_y, u_z)$ is the fluid velocity, $p$ is the pressure perturbation away from hydrostatic equilibrium, $T$ is the temperature perturbation away from radiative equilibrium (which has a temperature gradient $T_{0z}$ and an adiabatic temperature gradient $T_{ad,z}$), $\rho_0$ is the mean background density of the fluid in this region, $g = -g e_z$ is the gravity (and $e_z$ is the local vertical direction), and $\alpha$ is the thermal expansion coefficient. The quantities $T_{0z}$, $T_{ad,z}$, $\rho_0$, $g$, $\alpha$, $\nu$, and $\kappa_T$ are assumed to be constant.

Since the instability requires thermal diffusion to proceed, one of the major theoretical advances in this field has been the formal derivation of the so-called the Low Péclet Number (LPN) equations by Lignières (1999), who showed that whenever the local eddy-scale Péclet number is small, the evolution of the temperature perturbations is well approximated by solving

$$ u_z(T_{0z} - T_{ad,z}) = \kappa_T \nabla^2 T \quad (9) $$

for $T$ instead of evolving Equation (7). In this limit, the temperature perturbations and the vertical velocity are therefore slaved to one another, and $T$ can be solved when $u_z$ is known. This has a number of important consequences. For instance, it can easily be shown by taking a horizontal average of Equation (9) that temperature perturbations must have a zero mean if mass is conserved, which implies that the temperature stratification within the star remains unaffected by LPN turbulence. It is also obvious from Equation (9) that the system dynamics now only ever depend on the ratio $(T_{0z} - T_{ad,z})/\kappa_T$, which formally shows how a large thermal diffusivity can damp the effect of stratification. Finally, the total number of independent unknowns in the LPN limit is reduced by one compared with the full equations, which both saves computing time in DNSs and permits the derivation of interesting theoretical results, among which is the formal identification of criterion (2) as an energy stability criterion (Bischoff 2013; Garaud et al. 2015).

Using DNSs within the LPN limit, Prat & Lignières (2013, 2014), Garaud et al. (2015, 2017; Paper I hereafter), Prat et al. (2016), Garaud & Kulenthirarajah (2016) and Gagnier & Garaud (2018; Paper II hereafter) tested various aspects of Zahn's model. As reviewed in Papers I and II, the general conclusion of these studies is that Zahn's model performs remarkably well as long as its fundamental assumption, namely that the turbulence is local, applies. In particular, the works of Prat and collaborators agree with those of Garaud and collaborators in estimating the critical threshold for instability

in Equation (2) to be

$$ (J\mathrm{Pr})_c \approx 0.007. \quad (10) $$

Paper I also found that the turbulent mixing model given in Equation (5) fares relatively well whenever the local assumption is verified, with a coefficient

$$ C \approx 0.08 \quad (11) $$

for the turbulent diffusivity of a passive scalar and similarly for the turbulent viscosity. On the other hand, both Papers I and II demonstrated the need for an extension of Zahn's model in the limit where the local assumption is no longer valid (i.e., when the Zahn scale becomes commensurate with any system scale). Paper II also revealed that a region which does not satisfy the instability criterion (2) may still undergo substantial turbulent mixing if it is adjacent to one which does—a shear-driven analog to convective overshoot. Finally, Paper II demonstrated that in the strongly stratified limit, diffusive shear instabilities are intrinsically nonlinear and therefore exhibit hysteretic properties. As a result, a region that satisfies criterion (2) is not necessarily turbulent; whether it is or not depends on whether it is linearly unstable, and if not, then it depends on the evolutionary history of the shear layer.

In order to go beyond a strictly local description of turbulent transport in stratified shear flows, we propose in this paper a second-order turbulence closure model. This type of approach is well established, and first involves deriving exact evolution equations for the second-order moments of the equations (which include the Reynolds stress tensor and the turbulent buoyancy fluxes for instance). These equations in turn depend on third-order correlations, for which a closure model is then proposed. The closure model considered here is based on the one first put forward by Ogilvie (2003) and applied to unstratified shear-induced turbulence by Garaud & Ogilvie (2005) and to turbulent convection by Garaud et al. (2010). We now extend it to study stratified shear flows in the LPN limit. The construction of the model closely follows the steps outlined in Garaud et al. (2010), then proceeds to apply the LPN asymptotic expansion of Lignières (1999) to the closure equations. This is done in Section 2. Section 3 explores the local properties of the model and constrains its closure parameters against the existing DNSs of stratified plane Couette flow presented in Paper I. Section 4 compares the model predictions with the DNSs of body-forced stratified Kolmogorov flow presented in Garaud et al. (2015) and Garaud & Kulenthirarajah (2016). Section 5 compares them to the DNSs of Paper II and addresses model inadequacies through the introduction of an additional turbulent diffusivity term in the Reynolds stress equations. In Section 6, we propose a simple model for the turbulent eddy scale, which is a necessary ingredient of the closure model. We conclude in Section 7 by summarizing the findings of this paper and by discussing caveats and possible future extensions of the model.

## 2. A Second-order Turbulence Closure Model for Diffusive Shear Instabilities

### 2.1. Reynolds Stress Evolution Equations

We start from the momentum, thermal energy, and mass conservation Equations (6)–(8), now written for simplicity using Einstein's convention of summation over repeated
indices and generically nondimensionalized assuming a given velocity scale $U_0$, length scale $L_0$, pressure scale $p_0 = \rho_0 U_0^2$, and temperature scale $T_0 = (T_{0z} - T_{ad,z})L_0$. We have

$$ \partial_t u_i + \partial_k(u_i u_k) = -\partial_i p + \mathrm{RiT}\delta_{iz} + \frac{1}{\mathrm{Re}}\partial_{kk}^2 u_i (+\hat{F}(z)\delta_{ix}), \quad (12) $$

$$ \partial_t T + \partial_k(u_k T) + u_z = \frac{1}{\mathrm{Pe}}\partial_{kk}^2 T, \quad (13) $$

$$ \partial_k u_k = 0, \quad (14) $$

where the indices i and k can be x, y, or z (and similarly for the index j later on), and where

$$ \mathrm{Ri} = \frac{N^2 L_0^2}{U_0^2}, \quad \text{where } N^2 = \alpha g (T_{0z} - T_{ad,z}), \quad (15) $$

$$ \mathrm{Re} = \frac{U_0 L_0}{\nu}, \quad (16) $$

$$ \mathrm{Pe} = \frac{U_0 L_0}{\kappa_T}. \quad (17) $$

The quantities Re, Ri, and Pe constructed in this manner can be thought of, respectively, as the Richardson number, Reynolds number, and Péclet number for the flow. Since we will be comparing the closure model to different sets of DNSs that each use different normalizations, we leave the actual selection of $U_0$ and $L_0$ unspecified for now. Note that we have allowed for the existence of a nondimensional body force $\hat{F}(z)$ applied in the $e_x$ direction to drive the shear flow in Equation (12), but the latter is not required for wall-bounded shear flows (e.g., Paper I) or homogeneous shear flows in the shearing sheet approximation (Prat et al. 2016).

As in Ogilvie (2003), we now consider a suitable averaging procedure, which could be an average over realizations, or over a small length scale and/or short timescales, and which must commute with all differential operators ($\partial_i$ and $\partial_t$). This average is denoted with an overbar. Applying it to Equations (12)–(14) yields the mean equations:

$$ \partial_t \bar{u}_i + \partial_k(\overline{u_i u_k}) = -\partial_i \bar{p} + \mathrm{RiT}\delta_{iz} + \frac{1}{\mathrm{Re}}\partial_{kk}^2 \bar{u}_i (+\hat{F}(z)\delta_{ix}), \quad (18) $$

$$ \partial_t \bar{T} + \partial_k(\overline{u_k T}) + \bar{u}_z = \frac{1}{\mathrm{Pe}}\partial_{kk}^2 \bar{T}, \quad (19) $$

$$ \partial_k \bar{u}_k = 0. \quad (20) $$

We then separate all variables into mean and fluctuating parts such as $f = \bar{f} + f'$, where by construction $\bar{f}' = 0$. Using this decomposition in Equations (12)–(14) and subtracting the mean Equations (18)–(20), we obtain the evolution equations for the fluctuations:

$$ \begin{aligned} \partial_t u_i' + \bar{u}_k \partial_k u_i' + u_k' \partial_k \bar{u}_i + \partial_k(u_i' u_k') - \partial_k(\overline{u_i' u_k'}) \\ = -\partial_i p' + \mathrm{RiT}\delta_{iz} + \frac{1}{\mathrm{Re}}\partial_{kk}^2 u_i', \end{aligned} \quad (21) $$

$$ \begin{aligned} \partial_t T' + \bar{u}_k \partial_k T' + u_k' \partial_k \bar{T} + u_z' + \partial_k(u_k' T') - \partial_k(\overline{u_k' T'}) \\ = \frac{1}{\mathrm{Pe}}\partial_{kk}^2 T', \end{aligned} \quad (22) $$

$$ \partial_k u_k' = 0. \quad (23) $$

From these we can derive the exact equations for the evolution of the averages of quadratic quantities such as $R_{ij} = u_i' u_j'$ (the Reynolds stress), $F_i = u_i' T'$ (the turbulent temperature flux), and $Q = T'^2$ as follows:

$$ \begin{aligned} \partial_t \bar{R}_{ij} + \bar{u}_k \partial_k \bar{R}_{ij} + \bar{R}_{jk} \partial_k \bar{u}_j + \bar{R}_{ik} \partial_k \bar{u}_j - \mathrm{Ri} \bar{F}_i \delta_{iz} - \mathrm{Ri} \bar{F}_j \delta_{iz} \\ - \frac{1}{\mathrm{Re}} \partial_{kk}^2 \bar{R}_{ij} = -\overline{u_j' \partial_i p'} - \overline{u_i' \partial_j p'} - \overline{u_k' \partial_k R_{ij}} - \frac{2}{\mathrm{Re}} \overline{\partial_k u_i' \partial_k u_j'}, \end{aligned} \quad (24) $$

$$ \begin{aligned} \partial_t \bar{F}_i + \bar{u}_k \partial_k \bar{F}_i + \bar{F}_k \partial_k \bar{u}_i + \bar{R}_{ik} \partial_k \bar{T} \\ - \mathrm{Ri} \bar{Q} \delta_{iz} + \bar{R}_{iz} - \frac{1}{2} \left( \frac{1}{\mathrm{Re}} + \frac{1}{\mathrm{Pe}} \right) \partial_{kk}^2 \bar{F}_i \\ = - \overline{T' \partial_i p'} - \overline{T' \partial_k R_{ik} + u_i' \partial_k F_k} + \frac{1}{2} \left( \frac{1}{\mathrm{Re}} - \frac{1}{\mathrm{Pe}} \right) \\ \times (\partial_k(\overline{T' \partial_k u_i' - u_i' \partial_k T'})) - \left( \frac{1}{\mathrm{Re}} + \frac{1}{\mathrm{Pe}} \right) \overline{(\partial_k u_i' \partial_k T')}, \end{aligned} \quad (25) $$

$$ \begin{aligned} \partial_t \bar{Q} + \bar{u}_k \partial_k \bar{Q} + 2 \bar{F}_k \partial_k \bar{T} + 2 \bar{F}_z - \frac{1}{\mathrm{Pe}} \partial_{kk}^2 \bar{Q} \\ = - \overline{u_k' \partial_k \bar{Q}} - \frac{2}{\mathrm{Pe}} \overline{(\partial_k T')^2}. \end{aligned} \quad (26) $$

These equations are identical to those of Garaud et al. (2010) aside from the contribution associated with the assumed constant background temperature gradient which manifests itself in the terms $\bar{R}_{iz}$ in Equation (25) and $2\bar{F}_z$ in Equation (26). The left-hand side contains quantities that can be expressed exactly in terms of $\bar{R}_{ij}$, $\bar{F}_i$, $\bar{u}_i$, $\bar{T}$, and $\bar{Q}$ or their gradients. On the right-hand side are correlations involving the pressure fluctuations, triple correlations between fluctuating quantities, and dissipative terms involving microscopic diffusion. These terms cannot directly be expressed as functions of $\bar{R}_{ij}$, $\bar{F}_i$, and $\bar{Q}$ and must be modeled as part of the closure problem.

## 2.2. Proposed Second-order Closure Model

We adopt the simple closure model of Garaud et al. (2010; see also Ogilvie 2003; Garaud & Ogilvie 2005), applied to the set of Equations (24)–(26):

$$ \begin{aligned} \partial_t \bar{R}_{ij} + \bar{u}_k \partial_k \bar{R}_{ij} + \bar{R}_{ik} \partial_k \bar{u}_j + \bar{R}_{jk} \partial_k \bar{u}_i - \mathrm{Ri} \bar{F}_i \delta_{iz} - \mathrm{Ri} \bar{F}_j \delta_{iz} \\ - \frac{1}{\mathrm{Re}} \partial_{kk}^2 \bar{R}_{ij} = - \frac{C_1 \bar{R}^{1/2}}{L} \bar{R}_{ij} - \frac{C_2 \bar{R}^{1/2}}{L} \left( \bar{R}_{ij} - \frac{1}{3} \bar{R} \delta_{ij} \right) \\ - \frac{C_\nu \bar{R}_{ij}}{\mathrm{Re} L^2}, \end{aligned} \quad (27) $$

$$ \begin{aligned} \partial_t \bar{F}_i + \bar{u}_k \partial_k \bar{F}_i + \bar{F}_k \partial_k \bar{u}_i + \bar{R}_{ik} \partial_k \bar{T} - \mathrm{Ri} \bar{Q} \delta_{iz} + \bar{R}_{iz} \\ - \frac{1}{2} \left( \frac{1}{\mathrm{Re}} + \frac{1}{\mathrm{Pe}} \right) \partial_{kk}^2 \bar{F}_i = - \frac{C_6 \bar{R}^{1/2}}{L} \bar{F}_i \\ - \frac{1}{2L^2} \left( \frac{1}{\mathrm{Re}} + \frac{1}{\mathrm{Pe}} \right) C_{\nu\kappa} \bar{F}_i, \end{aligned} \quad (28) $$

$$ \begin{aligned} \partial_t \bar{Q} + \bar{u}_k \partial_k \bar{Q} + 2 \bar{F}_k \partial_k \bar{T} + 2 \bar{F}_z - \frac{1}{\mathrm{Pe}} \partial_{kk}^2 \bar{Q} \\ = - \frac{C_7 \bar{R}^{1/2}}{L} \bar{Q} - \frac{C_\kappa}{\mathrm{Pe} L^2} \bar{Q}, \end{aligned} \quad (29) $$
where $C_1, C_2, C_6, C_7, C_ν, C_k$, and $C_{νk}$ are assumed to be positive dimensionless universal constants, where $\bar{R}$ is the trace of the Reynolds stress tensor (which is also an estimate of the square of the turbulent rms velocity),

$$ \bar{R} = \bar{R}_{xx} + \bar{R}_{yy} + \bar{R}_{zz}, \qquad (30) $$

and $L$ is a characteristic length scale of the dominant turbulent eddies. The terms containing $C_1$ represent the decay of the turbulent kinetic energy through the turbulent cascade on the eddy turnover timescale of $L/\bar{R}^{1/2}$. The terms containing $C_2$ redistribute the turbulent energy through the different components of $\bar{R}_{ij}$ and model the expected isotropization of the turbulence on a similar timescale. The terms with $C_6$ and $C_7$ represent the turbulent decay of temperature fluctuations. The term containing $C_ν$ is introduced to model the viscous dissipation when an efficient turbulent cascade does not form at low Reynolds number, in which case the dissipation rate is directly proportional to the viscosity divided by the eddy scale $L$ squared. In a similar way, we introduce the last terms of Equations (25) and (26) containing $C_κ$ and $C_{νκ}$.

### 2.3. Closure Model in the LPN Limit

The equations derived so far could in principle be used to model shear-induced turbulence in high Péclet number stratified shear flows, but we do not recommend it. Indeed, stratified turbulence in that limit can be strongly influenced by gravity waves, whose long-range effects cannot be modeled within the context of this semi-local closure approach. However, gravity waves are not supported in the LPN limit (Lignières 1999), giving us hope that the closure model can be successfully applied here.

In order to derive a formal asymptotic limit of the closure model equations at low Péclet numbers, we apply the procedure of Lignières (1999), assuming a decomposition of the form $\bar{R}_{ij} = \bar{R}_{ij}^{(0)} + \text{Pe}\bar{R}_{ij}^{(1)} + O(\text{Pe}^2)$, $\bar{F}_i = \bar{F}_i^{(0)} + \text{Pe}\bar{F}_i^{(1)} + O(\text{Pe}^2)$, and $\bar{Q} = \bar{Q}^{(0)} + \text{Pe}\bar{Q}^{(1)} + O(\text{Pe}^2)$. To the lowest order in Pe, we have

$$ \partial_{kk}^2 \bar{F}_i^{(0)} = \frac{C_{\nu\kappa}}{L^2} \bar{F}_i^{(0)}, \quad \partial_{kk}^2 \bar{Q}^{(0)} = \frac{C_{\kappa}}{L^2} \bar{Q}^{(0)}, \qquad (31) $$

which do not have any bounded or periodic solutions. As a result of our particular model setup, $\bar{F}_i^{(0)} = \bar{Q}^{(0)} = 0$ for $i = x, y, z$. To the next order in Pe we have

$$ \partial_{kk}^2 \bar{Q}^{(1)} = \frac{C_\kappa}{L^2} \bar{Q}^{(1)}, \qquad (32) $$

which implies that $\bar{Q}^{(1)} = 0$ for the same reasons. This leaves

$$ \partial_{kk}^2 \bar{F}_i^{(1)} = 2\bar{R}_{iz}^{(0)} + \frac{C_{\nu\kappa}}{L^2} \bar{F}_i^{(1)}. \qquad (33) $$

Combining these equations with Equations (27)–(29) yields

$$ \begin{aligned} & \partial_t \bar{R}_{ij} + \bar{u}_k \partial_k \bar{R}_{ij} + \bar{R}_{ik} \partial_k \bar{u}_j + \bar{R}_{jk} \partial_k \bar{u}_i \\ & \quad - \text{RiPe}(\bar{f}_i \delta_{jz} + \bar{f}_j \delta_{iz}) - \frac{1}{\text{Re}} \partial_{kk}^2 \bar{R}_{ij} \\ & = -\frac{C_1 \bar{R}^{1/2}}{L} \bar{R}_{ij} - \frac{C_2 \bar{R}^{1/2}}{L} \left( \bar{R}_{ij} - \frac{1}{3} \bar{R} \delta_{ij} \right) - \frac{C_\nu \bar{R}_{ij}}{\text{Re} L^2}, \end{aligned} \qquad (34) $$

$$ \bar{R}_{iz} - \frac{1}{2} \partial_{kk}^2 \bar{f}_i = - \frac{1}{2L^2} C_{\nu\kappa} \bar{f}_i, \qquad (35) $$

where we have defined $f_i = F_i/\text{Pe}$ for simplicity. We recover the fundamental property of LPN flows, namely that their dynamics only depend on the product of the Richardson and Péclet numbers, RiPe.

In all that follows, we now take the average to imply a horizontal average. With this, we have $\bar{u}_z = 0$ by mass conservation, and assume for simplicity that $\bar{u}_y = 0$. Also, all of the horizontal derivatives disappear so

$$ \begin{aligned} & \partial_t \bar{R}_{xx} + 2\bar{R}_{xz} \partial_z \bar{u}_x - \frac{1}{\text{Re}} \partial_{zz}^2 \bar{R}_{xx} \\ & = -(C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{xx} + \frac{C_2 \bar{R}^{3/2}}{3L} - \frac{C_\nu \bar{R}_{xx}}{\text{Re} L^2}, \end{aligned} \qquad (36) $$

$$ \begin{aligned} & \partial_t \bar{R}_{yy} - \frac{1}{\text{Re}} \partial_{zz}^2 \bar{R}_{yy} = -(C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{yy} \\ & \quad + \frac{C_2 \bar{R}^{3/2}}{3L} - \frac{C_\nu \bar{R}_{yy}}{\text{Re} L^2}, \end{aligned} \qquad (37) $$

$$ \begin{aligned} & \partial_t \bar{R}_{zz} - 2\text{RiPe}\bar{f}_z - \frac{1}{\text{Re}} \partial_{zz}^2 \bar{R}_{zz} \\ & = -(C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{zz} + \frac{C_2 \bar{R}^{3/2}}{3L} - \frac{C_\nu \bar{R}_{zz}}{\text{Re} L^2}, \end{aligned} \qquad (38) $$

$$ \begin{aligned} & \partial_t \bar{R}_{xz} + \bar{R}_{zz} \partial_z \bar{u}_x - \text{RiPe}\bar{f}_x - \frac{1}{\text{Re}} \partial_{zz}^2 \bar{R}_{xz} \\ & = -(C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{xz} - \frac{C_\nu \bar{R}_{xz}}{\text{Re} L^2}, \end{aligned} \qquad (39) $$

$$ \begin{aligned} & \partial_t \bar{R}_{yz} - \text{RiPe}\bar{f}_y - \frac{1}{\text{Re}} \partial_{zz}^2 \bar{R}_{yz} = -(C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{yz} \\ & \quad - \frac{C_\nu \bar{R}_{yz}}{\text{Re} L^2}, \end{aligned} \qquad (40) $$

$$ \begin{aligned} & \partial_t \bar{R}_{xy} + \bar{R}_{yz} \partial_z \bar{u}_x - \frac{1}{\text{Re}} \partial_{zz}^2 \bar{R}_{xy} = -(C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{xy} \\ & \quad - \frac{C_\nu \bar{R}_{xy}}{\text{Re} L^2}, \end{aligned} \qquad (41) $$

$$ \bar{R}_{xz} - \frac{1}{2} \partial_{zz}^2 \bar{f}_x = - \frac{1}{2L^2} C_{\nu\kappa} \bar{f}_x, \qquad (42) $$

$$ \bar{R}_{yz} - \frac{1}{2} \partial_{zz}^2 \bar{f}_y = - \frac{1}{2L^2} C_{\nu\kappa} \bar{f}_y, \qquad (43) $$

$$ \bar{R}_{zz} - \frac{1}{2} \partial_{zz}^2 \bar{f}_z = - \frac{1}{2L^2} C_{\nu\kappa} \bar{f}_z. \qquad (44) $$

Note how the temperature correlations $\bar{Q}$ and the original model constants $C_6, C_7$, and $C_\kappa$ drop out of the LPN closure model. Also note that the $\bar{f}_y$ and $\bar{R}_{yz}$ equations decouple from the others. Since they are not forced, these quantities must eventually always decay away. In what follows, we now assume them to be zero. In addition, the quantity $\bar{R}_{xy}$ is not involved in any equation aside from its own evolution equation. We can therefore ignore it until needed for some particular purpose. The remaining coupled system of equations reduces to Equations (36)–(39) as well as to Equations (42) and (44), together with the mean flow Equation (18), where $\bar{T} = 0$ in the LPN limit (see Section 1).
### 3. Local Properties of the Closure Model, and Parameter Fitting

#### 3.1. Local Properties of the Model

It is informative to investigate what the closure model predictions are for homogeneous stratified shear flows. To do so, we assume that the background shear $\partial_z \bar{u}_x \equiv S$ is known and constant, and ignore all other derivatives in the closure equations accordingly. In this section, we ignore the body force and select the following nondimensionalization: $S^{-1}$ as unit time, $L_z$ as unit length (where $L_z$ is the height of the computational domain considered), and $L_z(T_{0z} - T_{0z}^{\text{ad}})$ as unit temperature. This redefines the Richardson, Péclet, and Reynolds numbers to be $R_i \equiv J = N^2/S^2$, $Pe \equiv Pe_S = SL_z^2/\kappa_T$, and $Re \equiv Re_S = SL_z^2/\nu$. The closure equations then reduce to

$$2\bar{R}_{xz}s = -\frac{C_1 + C_2}{L}\bar{R}^{1/2}\bar{R}_{xx} + \frac{C_2}{3L}\bar{R}^{3/2} - \frac{1}{\text{Re}_S}\frac{C_\nu}{L^2}\bar{R}_{xx}, \quad (45)$$

$$0 = -\frac{C_1 + C_2}{L}\bar{R}^{1/2}\bar{R}_{yy} + \frac{C_2}{3L}\bar{R}^{3/2} - \frac{1}{\text{Re}_S}\frac{C_\nu}{L^2}\bar{R}_{yy}, \quad (46)$$

$$-2J\text{Pe}_S\bar{f}_z = -\frac{C_1 + C_2}{L}\bar{R}^{1/2}\bar{R}_{zz} + \frac{C_2}{3L}\bar{R}^{3/2} - \frac{1}{\text{Re}_S}\frac{C_\nu}{L^2}\bar{R}_{zz}, \quad (47)$$

$$\bar{R}_{zz}s - J\text{Pe}_S\bar{f}_x = -\frac{C_1 + C_2}{L}\bar{R}^{1/2}\bar{R}_{xz} - \frac{1}{\text{Re}_S}\frac{C_\nu}{L^2}\bar{R}_{xz}, \quad (48)$$

$$\bar{R}_{xz} = -\frac{C_{\nu\kappa}}{2L^2\bar{f}_x}, \quad (49)$$

$$\bar{R}_{zz} = -\frac{C_{\nu\kappa}}{2L^2\bar{f}_z}, \quad (50)$$

where $L$ is now implicitly given in units of $L_z$, and $s = \text{sign}(S)$. The two flux equations can straightforwardly be solved for $\bar{f}_x$ and $\bar{f}_z$, and substituted into the Reynolds stress equations. Furthermore, we can construct an equation for $\bar{R}$ by summing the $\bar{R}_{xx}$, $\bar{R}_{yy}$, and $\bar{R}_{zz}$ equations. Once this is done, both $\bar{R}_{xx}$ and $\bar{R}_{yy}$ decouple from the system, which can be reduced to three coupled equations only (for $\bar{R}$, $\bar{R}_{zz}$, and $\bar{R}_{xz}$):

$$4J\text{Pe}_S \frac{L^2}{C_{\nu\kappa}} \bar{R}_{zz} = -\frac{C_1 + C_2}{L} \bar{R}^{1/2} \bar{R}_{zz} + \frac{C_2}{3L} \bar{R}^{3/2} - \frac{1}{\text{Re}_S} \frac{C_\nu}{L^2} \bar{R}_{zz}, \quad (51)$$

$$2\bar{R}_{xz}s + 4J\text{Pe}_S \frac{L^2}{C_{\nu\kappa}}\bar{R}_{zz} = -\frac{C_1}{L}\bar{R}^{3/2} - \frac{1}{\text{Re}_S}\frac{C_\nu}{L^2}\bar{R}, \quad (52)$$

$$\bar{R}_{zz}s + 2J\text{Pe}_S \frac{L^2}{C_{\nu\kappa}}\bar{R}_{xz} = -\frac{C_1 + C_2}{L}\bar{R}^{1/2}\bar{R}_{xz} - \frac{1}{\text{Re}_S}\frac{C_\nu}{L^2}\bar{R}_{xz}. \quad (53)$$

We see that two parameters naturally emerge, namely $\sigma = J\text{Pe}_S L^2/C_{\nu\kappa}$, which controls the stratification, and $\epsilon = C_\nu/\text{Re}_S L^2$, which controls macroscopic viscous damping on the scale $L$. Furthermore, if we define $\lambda = \bar{R}^{1/2}/L$, $r_{zz} = \bar{R}_{zz}/\bar{R}$, and similarly for $r_{xz}$, we get the simpler-looking equations

$$4\sigma r_{zz} = -(C_1 + C_2)\lambda r_{zz} + \frac{C_2}{3}\lambda - \epsilon r_{zz}, \quad (54)$$

$$2r_{xz}s + 4\sigma r_{zz} = -C_1\lambda - \epsilon, \quad (55)$$

$$r_{zz}s + 2\sigma r_{xz} = -(C_1 + C_2)\lambda r_{xz} - \epsilon r_{xz}, \quad (56)$$

which can be combined into a single cubic equation for $\lambda$:

$$\begin{aligned} & \frac{2C_2}{3}\lambda \\ = & \left\{ 4\sigma \left[ \left(C_1 + \frac{C_2}{3}\right)\lambda + \epsilon \right] + (C_1\lambda + \epsilon)((C_1 + C_2)\lambda + \epsilon) \right\} \\ & \times [2\sigma + (C_1 + C_2)\lambda + \epsilon]. \end{aligned} \quad (57)$$

For general values of the parameters $\sigma$ and $\epsilon$, this equation must be solved numerically, and from there, the values of $\bar{R}$ and of all the components of the stress tensor can be recovered. In some limits, analytical solutions also exist.

For instance, the unstratified limit ($\sigma=0$) was studied by Garaud & Ogilvie (2005). In that case,

$$\frac{2C_2}{3}\lambda = (C_1\lambda + \epsilon)((C_1 + C_2)\lambda + \epsilon)^2, \quad (58)$$

the solutions of which can be computed numerically for given values of $\epsilon$, $C_1$, and $C_2$ and have the following asymptotic approximation in the limit of small $\epsilon$ (or equivalently, when $\text{Re}_S \to \infty$):

$$\lambda_1 = \frac{3\epsilon^3}{2C_2} + O(\epsilon^4), \quad (59)$$

$$\lambda_2 = \left( \frac{2C_2}{3C_1} \frac{1}{(C_1 + C_2)^2} \right)^{1/2} - \frac{3C_1 + C_2}{2C_1(C_1 + C_2)}\epsilon + O(\epsilon^2), \quad (60)$$

$$\lambda_3 = -\left(\frac{2C_2}{3C_1}\frac{1}{(C_1 + C_2)^2}\right)^{1/2} - \frac{3C_1 + C_2}{2C_1(C_1 + C_2)}\epsilon + O(\epsilon^2). \quad (61)$$

The negative solution $\lambda_3$ can be discarded on the basis that $\lambda$ has to be positive. As shown by Garaud & Ogilvie (2005), $\lambda_1$ corresponds to a repelling stationary turbulent state (and can therefore also be discarded) while $\lambda_2$ corresponds to the attracting stationary turbulent state that we seek. The solutions are shown in Figure 1(a), for fiducial parameters $C_1 = 0.41$ and $C_2 = 0.54$.

In the case where both $\sigma = 0$ and $\epsilon = 0$, the stress tensor components corresponding to the root $\lambda_2$ have simple exact analytical expressions (Garaud & Ogilvie 2005):

$$\begin{aligned} \lambda &= \frac{\bar{R}^{1/2}}{L} = \frac{1}{C_1 + C_2} \sqrt{\frac{2C_2}{3C_1}}, \\ r_{xz} &= -\frac{sC_1}{2(C_1 + C_2)} \sqrt{\frac{2C_2}{3C_1}}, \quad r_{xx} = \frac{3C_1 + C_2}{3(C_1 + C_2)}, \\ r_{yy} &= r_{zz} = \frac{C_2}{3(C_1 + C_2)}. \end{aligned} \quad (62)$$

Notably, we see that all of the ratios $r_{ij} = \bar{R}_{ij}/\bar{R}$ depend only on the ratio of the model constants $C_2/C_1$. This is an important property of the closure model that can both be tested numerically and used to constrain $C_2/C_1$ (see Section 3.2).

When $\sigma \neq 0$, Equation (57) can only be solved numerically. Three branches of solutions exist, and as before, only two of them have $\lambda > 0$ and only one corresponds to an attracting state. They are shown in Figure 1(a) for different values of $\sigma$, as
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left: numerical solutions for $\lambda_1$ (repelling turbulent branch, dashed line) and $\lambda_2$ (attracting turbulent branch, solid line) as a function of $\epsilon^{-1}$ (which is proportional to the Reynolds number), for fiducial parameters $C_1 = 0.41$ and $C_2 = 0.54$ (see Section 3.2), and for various values of the stratification parameter $\sigma$ (the value of $\sigma$ is marked near its corresponding curve). Right: plot of $\epsilon\sigma_c(\epsilon)$ as a function of $\epsilon^{-1}$ for $C_1 = 0.41$ and $C_2 = 0.54$. This quantity is related to the critical value of JPr in diffusive stratified shear flows (see the text for details).</figcaption></figure>

a function of $\epsilon^{-1}$ (which is proportional to the Reynolds number of the flow). We see that the transition to turbulence occurs through a finite amplitude instability rather than a linear instability at $\sigma = 0$ (Garaud & Ogilvie 2005), and this remains the case as $\sigma$ increases. We also see that the critical Reynolds number above which turbulent solutions exist increases with $\sigma$. At fixed values of $\epsilon$, on the other hand, turbulent solutions disappear when $\sigma$ exceeds a critical value $\sigma_c(\epsilon)$, a result that is consistent with the notion that a shear flow can be stabilized by sufficiently large stratification even in the LPN limit.

In fact, it is interesting to note that

$$ \epsilon\sigma = \frac{C_{\nu}}{C_{\nu\kappa}} \text{JPr}, \qquad (63) $$

so this product is important for several reasons. First, it no longer depends on the (so far arbitrary) selection of the eddy scale $L$ (on which both $\sigma$ and $\epsilon$ depend), and second, the quantity JPr was found to be the relevant bifurcation parameter for stratified shear flows in the LPN limit, with flows having $JPr < (JPr)_c \simeq 0.007$ being unstable and flows with $JPr > (JPr)_c$ being stable to diffusive shear instabilities (see Section 1). For this reason, we show in Figure 1(b) the product $\epsilon\sigma_c(\epsilon)$ computed from the numerical solutions shown in Figure 1(a). We see that this product is not constant but, instead, exhibits a well-defined maximum at intermediate values of $\epsilon^{-1}$ and tends to zero in the inviscid limit. The fact that $\epsilon\sigma_c(\epsilon) \to 0$ in the limit of small $\epsilon^{-1}$ is not actually surprising: a shear flow should indeed become stable even in the unstratified limit for sufficiently large viscosity. On the other hand, the fact that $\epsilon\sigma_c(\epsilon) \to 0$ in the limit of small viscosity (large $\epsilon^{-1}$) is somewhat more puzzling, since we would indeed expect that this quantity should tend to a constant in the inviscid limit, at least in the light Zahn model (Zahn 1974). As we shall demonstrate in Section 3.2, however, this apparent discrepancy is not actually problematic, and some elements of Zahn's model can indeed be recovered.

Finally, the strictly inviscid limit $\epsilon = 0$ for arbitrary values of $\sigma$ is easy to solve analytically, since we have

$$ \frac{2C_2}{3}\lambda = \lambda \left[ 4\sigma \left( C_1 + \frac{C_2}{3} \right) + C_1(C_1 + C_2)\lambda \right] \left[ 2\sigma + (C_1 + C_2)\lambda \right]. \quad (64) $$

One of the roots is $\lambda_1 = 0$, while the other two satisfy the quadratic

$$ \frac{2C_2}{3} = \left[ 4\sigma \left( C_1 + \frac{C_2}{3} \right) + C_1(C_1 + C_2)\lambda \right] [2\sigma + (C_1 + C_2)\lambda], \quad (65) $$

whose solutions are

$$ \lambda_{2,3} = - \left[ \left( 3 + \frac{2C_2}{3C_1} \right) \pm \sqrt{ \left( 1 + \frac{2C_2}{3C_1} \right)^2 + \frac{2C_2}{3C_1} \frac{1}{\sigma^2} } \right] \frac{\sigma}{C_1 + C_2}. \quad (66) $$

While this always has real solutions, a necessary condition for at least one of these solutions to be positive is

$$ 4 \left( \frac{3C_1}{C_2} + 1 \right) \sigma^2 \le 1. \quad (67) $$

This states that turbulence dominated by eddies of scale $L$ can only be sustained in inviscid LPN stratified shear flows provided $\sigma < \sigma_c(\epsilon = 0) = 0.5(3C_1/C_2 + 1)^{-1/2} \simeq 0.276$ for the fiducial parameters $C_1 = 0.41$ and $C_2 = 0.54$ (see Section 3.2). We see in Figure 1(a) (see also Figure 6) that $\lambda \to 0$ as $\sigma \to \sigma_c$, and that the critical value of $\epsilon^{-1}$ above which turbulent solutions exist tends to infinity in the same limit. In other words, turbulent solutions become weaker as $\sigma \to \sigma_c$ and more difficult to sustain. Above $\sigma_c$, they disappear entirely. Interestingly, if we define the Péclet number of eddies on scale $L$ as $Pe_L = L^2Pe_s$, then we can express the condition
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Numerical solutions for $\lambda(\sigma)$ in the limit of $\epsilon \to 0$ for our parameters $C_1 = 0.41$ and $C_2 = 0.54$ (solid line). Also shown is the curve $\sqrt{2C_2/3C_1}(1 - \sigma/\sigma_c)/(C_1 + C_2)$ (dashed line), which is a relatively good approximation to the true solution for these parameters.</figcaption></figure>

for instability in the inviscid limit as

$$JPe_L \le \frac{C_{\nu\kappa}}{2} \left( \frac{3C_1}{C_2} + 1 \right)^{-1/2} . \qquad (68)$$

This criterion is reminiscent of the findings of Townsend (1958) and Dudis (1974) for the linear instability of stratified shear flows in the optically thin regime, even though it is applied here to optically thick homogeneous stratified shear flows that are linearly stable (Knobloch 1984) but nonlinearly unstable. It is also at the heart of Zahn's general nonlinear criterion for diffusive shear instabilities (see Section 1 and Paper I for a review of the topic). It is therefore reassuring to see that our closure model captures this essential feature of diffusive shear instabilities.

The positive solution of Equation (65) is shown in Figure 2 for $C_1 = 0.41$, $C_2 = 0.54$, and $0 < \sigma < \sigma_c(\epsilon=0)$. It happens to be relatively well approximated by the linear function $\lambda \approx \sqrt{2C_2/3C_1}(1 - \sigma/\sigma_c)/(C_1 + C_2)$ for $0 \le \sigma \le \sigma_c$, and this formula can then be used to reconstruct all of the components of the stress tensor if an approximate solution to the local closure model is all that is needed. Note that this linear function was chosen because it satisfies $\lambda = 0$ when $\sigma = \sigma_c$, and $\lambda = \sqrt{2C_2/3C_1}/(C_1 + C_2)$ when $\sigma = 0$, which are exact solutions for these particular values of $\sigma$.

## 3.2. Comparison with Numerical Simulations of Plane Couette Flow and Parameters Fitting

The results of the local model can be used to constrain the closure parameters $C_1, C_2, C_\nu$, and $C_{\nu\kappa}$, and to gain better insight into possible strategies for selecting the eddy scale $L$. Ideally, we should constrain the parameters by comparing the model with simulations of homogeneous stratified shear flows, such as the ones produced by Prat & Lignières (2013, 2014) or Prat et al. (2016). However, these papers do not report the values of all components of the stress tensor measured in their simulations and therefore cannot be used for our purpose. Instead, we elect to compare the model with the results of the DNSs of stratified plane Couette flow of Paper I. While such a setup was not originally designed to produce homogeneous

stratified turbulence (because it is bounded by the top and bottom plates), we have found that many of the turbulence properties in the central region of the domain, far from the boundaries, are quantitatively consistent with those measured in the shearing box simulations of Prat et al. (2016; see Paper I for details). In what follows, we first briefly summarize the DNSs of Paper I and explain how various quantities of interest, such as the Reynolds stresses, are extracted and then used to constrain the closure parameters.

### 3.2.1. Description of the DNSs from Paper I

Paper I presents DNSs of the stratified plane Couette flow in the LPN limit. This setup is boundary driven: (1) the shear is driven between two no-slip plane-parallel plates moving at velocities $\pm \Delta U/2$, and (2) a background temperature gradient is maintained by holding the top and bottom plates at temperatures $\pm \Delta T/2$. The vertical distance between the two plates is $L_z$, and the domain is assumed to be periodic in the horizontal direction. The system dynamics are evolved using the nondimensional LPN equation

$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = - \nabla p + \mathrm{Ri}_C \mathrm{Pe}_C \nabla^{-2} u_z \mathbf{e}_z + \frac{1}{\mathrm{Re}_C} \nabla^2 \mathbf{u} \quad \text{with} \quad \nabla \cdot \mathbf{u} = 0, \qquad (69)$$

where $\mathrm{Ri}_C = \alpha g L_z \Delta T / \Delta U^2$, $\mathrm{Pe}_C = \Delta U L_z / \kappa_T$, and $\mathrm{Re}_C = \Delta U L_z / \nu$ are the global Richardson, Péclet, and Reynolds numbers for this problem. In all cases, the solution is evolved numerically in time until a statistically stationary state is reached, at which point the simulation is analyzed to extract the quantities of interest (see Paper I for details of the code). In that state, the shear in the bulk of the domain is nearly constant aside from two thin viscous boundary layers located in the immediate vicinity of each plate. Paper I already reported the value of the mean shear in the bulk of the domain (away from the boundaries) $S$ (reported as $\hat{S}$ and measured in units of $\Delta U/L_z$), the component of the Reynolds stress tensor $\bar{R}_{xz}$ (reported as $\langle uw \rangle_i$ and measured in units of $\Delta U^2$), as well as the typical eddy scale $l_e$ (in units of $L_z$) in Table 1. See Paper I for details on the method used to obtain these bulk quantities. We have reanalyzed the same simulations to extract in addition $\bar{R}_{xx}$, $\bar{R}_{yy}$, and $\bar{R}_{zz}$ in the bulk of the domain for the two sets of runs with the highest Reynolds numbers only ($\mathrm{Re}_C = 9 \times 10^4$ and $1.2 \times 10^5$).

Knowing these three components, we can then calculate $\bar{R} = \bar{R}_{xx} + \bar{R}_{yy} + \bar{R}_{zz}$ and then the corresponding ratios $r_{xz}, r_{xx}, r_{yy}$, and $r_{zz}$ (see Section 3.1). Finally, the value of $JPe_S$ corresponding to the shear in the middle of the domain is calculated using the formula $JPe_S = \mathrm{Ri}_C \mathrm{Pe}_C / S$.

### 3.2.2. Constraining the Model Parameters

We begin by comparing the results of DNSs in the very weakly stratified limit with the analytical formulas given in Equation (62). As discussed in Section 3.1, the local model predictions for the ratios $r_{xz}, r_{xx}, r_{yy}$, and $r_{zz}$ are independent of the eddy scale $L$ and only depend on the ratio $C_2/C_1$, which provides a way of testing the model and constraining that ratio using several independent measurements. Figure 3 shows $r_{xz}, r_{xx}, r_{yy}$, and $r_{zz}$ as a function of $JPe_S$ for the two highest
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Variation of the ratios $r_{xz}$, $r_{xx}$, $r_{yy}$, and $r_{zz}$ with $J\text{Pe}_S$ in the two highest Reynolds number sets of simulations presented in Paper I. The dashed lines correspond to the $\text{Re}_C = 9 \times 10^4$ simulations while the solid lines correspond to the $\text{Re}_C = 1.2 \times 10^5$ simulations. The error bars reflect the measured rms fluctuations of each of the $\bar{R}_{ij}$; only those for $\text{Re}_C = 1.2 \times 10^5$ are shown to avoid crowding the figure.</figcaption></figure>

Reynolds number sets of simulations available. We see that the results of the two sets are consistent with one another, confirming that the Reynolds number is indeed sufficiently high to be in the inviscid limit for these runs.

As predicted by the local closure model, $r_{yy}$ is indeed very close to $r_{zz}$ in the limit of weak stratification ($J\text{Pe}_S \ll 1$). We also see that in the same limit, $r_{xx} \simeq 0.62$, $r_{zz} \simeq 0.19$, and $r_{xz} \simeq 0.15$. Comparing these values with their analytical predictions would imply that $C_2/C_1$ should be roughly equal to 1.32, 1.32, and 1.34, respectively, three values that are remarkably consistent with one another. This is an encouraging indication that the closure model is behaving appropriately in this unstratified limit. In order to constrain $C_1$ and $C_2$ individually, we use the results of Garaud & Ogilvie (2005), who noted that the quantity

$$ \kappa = \frac{2}{C_1} \left[ \frac{C_1 C_2}{6(C_1 + C_2)^2} \right]^{3/4} \quad (70) $$

is the well-known von Kármán constant of the Prandtl–von Kármán universal velocity profile for wall-bounded shear flows. This constant has been measured in a number of laboratory experiments, notably those of Zagarola & Smits (1998), and found to be $\kappa = 0.436$. Using $C_2/C_1 = 1.32$ in this formula, we finally find that $C_1 \simeq 0.41$ and so $C_2 \simeq 0.54$. The value of $C_1$ found is remarkably close to that originally proposed by Garaud & Ogilvie (2005), while $C_2$ is smaller than in Garaud & Ogilvie (2005) by 10% only. In what follows, we therefore adopt $C_1 = 0.41$ and $C_2 = 0.54$ as our fiducial parameters.

One would also like to compare the model prediction for $\bar{R}$ (given in Equation (62)) with the available data. Unlike the ratios $r_{ij}$, however, $\bar{R}$ depends on the assumed eddy scale $L$ that has so far remained unspecified. For this reason, we use a different approach: instead of picking an arbitrary model for $L$ and comparing the predicted values of $\bar{R}$ with the data, we use Equation (62) to compute what value of $L$ would be required to reproduce the values of $\bar{R}$ measured in the DNSs for low stratification. In other words, we compute the quantity $L_{\text{fit}} = (C_1 + C_2)(3C_1/2C_2)^{1/2}\bar{R}^{1/2}$ from

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of the value of $L_{\text{fit}}$ required to fit the data from the DNSs of Paper I for $\bar{R}^{1/2}$ (blue), with the measured eddy scale $l_e$ (red). The dashed lines correspond to the $\text{Re}_C = 9 \times 10^4$ simulations while the solid lines correspond to the $\text{Re}_C = 1.2 \times 10^5$ simulations. The error bars on $L_{\text{fit}}$ reflect the measured rms fluctuations of $\bar{R}$; only those for $\text{Re}_C = 1.2 \times 10^5$ are shown to avoid crowding the figure. Error bars on $l_e$ are not easy to measure, but would be of the same order as those measured for $L_{\text{fit}}$. Note that the comparison between $l_e$ and $L_{\text{fit}}$ is only meaningful in the limit $J\text{Pe}_S \ll 1$, where Equation (62) applies. We find that in this limit $L_{\text{fit}} \simeq 1.5l_e$.</figcaption></figure>

the simulations, using the fiducial values of $C_1$ and $C_2$ selected earlier. The results are shown in Figure 4 for the two sets of DNSs from Paper I with the highest Reynolds numbers and compared with the actual eddy scale $l_e$ measured from the same DNSs (see Table 1 of Paper I). Note that the comparison is only meaningful in the limit where $J\text{Pe}_S \ll 1$, for which Equation (62) applies. We see that the value $L_{\text{fit}}$ required to fit the $\bar{R}$ data is a little larger than $l_e$ by a factor of 1.5 in this limit. The discrepancy is not surprising, given that the meanings of $L$ and $l_e$ are very different: one is created from an energy cascade timescale in a closure model, and the other is measured from DNSs using an autocorrelation length scale (see Paper I). As such, we do not expect them to be identical a priori, but we do expect them to remain proportional to one another as parameters are varied. In the remainder of this section, and for the purpose of calibrating the remaining model parameters $C_\nu$ and $C_{\kappa\nu}$, we now take $L = 1.5l_e$, where $l_e$ is the measured eddy scale from DNSs. This assumption is successfully verified in Section 4.

In order to constrain $C_{\kappa\nu}$, we must now look at the limit where $\sigma$ is not asymptotically small. Assuming as before that viscosity is negligible, we have

$$ \begin{aligned} \frac{\bar{R}^{1/2}}{L} &= \lambda(\sigma) \\ &= -\left[ \left( 3 + \frac{2C_2}{3C_1} \right) - \sqrt{ \left( 1 + \frac{2C_2}{3C_1} \right)^2 + \frac{2C_2}{3C_1} \frac{1}{\sigma^2} } \right] \frac{\sigma}{C_1 + C_2}, \end{aligned} \quad (71) $$

where we have selected the solution that is not always negative. The right-hand side of this expression is the function $\lambda(\sigma)$ shown in Figure 2, which only depends on $\sigma = J\text{Pe}_S L^2/C_{\kappa\nu}$, where $J\text{Pe}_S L^2$ is known from the data (assuming that $L = 1.5l_e$, where $l_e$ is measured in the DNSs). This implies that it contains only one unknown, namely $C_{\kappa\nu}$. Meanwhile, the left-hand side can be directly measured from the available data too. We can then vary $C_{\kappa\nu}$ (which varies $\sigma$) to find the best fit between the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of the left-hand side and right-hand side of Equation (71), using the data from the DNSs of Paper I to compute $\bar{R}^{1/2}/L$, with $L$ taken to be $1.5l_e$. The dashed lines correspond to the $Re_c = 9 \times 10^4$ simulations while the solid lines correspond to the $Re_c = 1.2 \times 10^5$ simulations. The error bars on $\bar{R}^{1/2}/L$ reflect the measured rms fluctuations of $\bar{R}$; only those for $Re_c = 1.2 \times 10^5$ are shown to avoid crowding the figure. The best fit, shown here, is obtained with $C_{\nu\kappa} = 10$.</figcaption></figure>

left- and right-hand sides of Equation (71). The results are shown in Figure 5, with a best fit obtained for $C_{\nu\kappa} \simeq 10$. This is fairly significantly larger than the value of $C_{\nu\kappa} \simeq 6$ estimated by Garaud et al. (2010) by fitting their related closure model against DNSs of Rayleigh-Bénard convection but is still within a factor of 2 of the latter.

The most difficult parameter to constrain is $C_\nu$, since the latter is only relevant when viscosity becomes relevant, a limit we have so far avoided. Comparing the model to the lowest Reynolds number DNSs of Paper I is not a good idea, because the influence of the viscous boundary layers (where the local model breaks down) is larger for lower Reynolds numbers. Instead, we calibrate $C_\nu$ (which appears in the definition of $\epsilon$) to ensure that the closure model correctly predicts the critical value of the stratification above which the shear flow becomes stable in the DNSs.

As explained by Zahn (1974) and seen in Figure 4, the dominant eddy scale of diffusive shear instabilities decreases as the stratification increases (see Paper I for more on this topic). For sufficiently large stratification, that scale becomes viscously damped, at which point the shear is completely stabilized. In Paper I, we showed that this happens in DNSs when $JPr = (JPr)_c \simeq 0.007$, while from the closure model properties we find that this happens when $\sigma \sim \sigma_c(\epsilon)$. Earlier, we also showed that $\epsilon\sigma = JPrC_\nu/C_{\nu\kappa}$. As a result, we can calibrate $C_\nu$ to ensure that the numerically determined stability limit $JPr \simeq 0.007$ is consistent with the closure model limit $\sigma = \sigma_c(\epsilon)$ (or equivalently, $\epsilon\sigma = JPrC_\nu/C_{\nu\kappa} = \epsilon\sigma_c(\epsilon)$). This is done in Figure 6, which shows $JPr$ as a function of $\sigma$ and compares several curves: the numerically determined threshold ($JPr)_c$ from Paper I, the actual data points from the two sets of DNSs with the highest Reynolds numbers of Paper I, and the critical threshold for stability in our closure model, namely $\sigma\epsilon_c(\sigma) C_{\nu\kappa}/C_\nu$, for $C_{\nu\kappa} = 10$ and different values of $C_\nu$. The value of $C_\nu$ that appropriately captures the correct threshold from the numerical simulations is $C_\nu \sim 15$. This value is very much consistent with the one originally estimated by Garaud & Ogilvie (2005), namely $C_\nu \sim 12.5$, which is remarkable given

the fact that the latter was determined using a completely different approach (by fitting the velocity profile for wall-bounded shear flows).

To conclude, in what follows we therefore adopt the parameters $C_1 = 0.41$, $C_2 = 0.54$, $C_{\nu\kappa} = 10$, and $C_\nu = 15$. Note that we have refrained from providing error bars on the estimated values for the model parameters, because these error bars would only be meaningful in the context of the specific method we used to fit them. In other words, had we selected a different method or a different data set, we would likely have obtained somewhat different values of the parameters, with their own error bars. The only meaningful way to estimate the errors on $C_1$, $C_2$, $C_{\nu\kappa}$, and $C_\nu$ is to fit the model to many different types of systems and to see how different the estimated parameters are. We have already seen that the values of $C_1$, $C_2$, and $C_\nu$ fitted against the Paper I DNSs are within ~20% of the corresponding values measured by Garaud & Ogilvie (2005). Meanwhile, $C_{\nu\kappa}$ is within ~50% of the value proposed by Garaud et al. (2010). Hence, 20% (respectively 50%) could be a good estimate of the size of the error on the model parameters $C_1$, $C_2$, and $C_\nu$ (respectively $C_{\nu\kappa}$).

## 4. Comparison of the Closure Model Predictions with DNSs of Stratified Kolmogorov Flows in a Statistically Stationary State

Having estimated the closure model parameters $C_1$, $C_2$, $C_\nu$, and $C_{\nu\kappa}$ from DNSs of stratified plane Couette flow, we now compare the predictions of the model with results from DNSs of stratified Kolmogorov flows (i.e., spatially sinusoidal flows) of Garaud & Kulenthirarajah (2016). These sets of simulations solve the LPN equations

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = - \nabla p + \mathrm{Ri}_F \mathrm{Pe}_F \nabla^{-2} u_z e_z + \frac{1}{\mathrm{Re}_F} \nabla^2 \mathbf{u} + \hat{F}(z) e_x \text{ with } \nabla \cdot \mathbf{u} = 0, \qquad (72) $$

where the nondimensionalization used is based on the amplitude $F_0$ and wavenumber $k$ of the sinusoidal forcing function (see Garaud & Kulenthirarajah 2016 for details). The unit length $L_0 = k^{-1}$ and velocity $U_0 = U_F = (F_0/k\rho_0)^{1/2}$ are used to create appropriate Richardson, Péclet, and Reynolds numbers as $\mathrm{Ri}_F = N^2/k^2 U_F^2$, $\mathrm{Pe}_F = U_F/k\kappa_T$, and $\mathrm{Re}_F = U_F/k\nu$. The non-dimensional forcing becomes $\hat{F}(z) = \sin(z)$, and the nondimensional height of the domain is $2\pi$. The perturbations are triply periodic. All runs available are summarized in Table 3 of Garaud & Kulenthirarajah (2016). In what follows, we focus once again on the highest Reynolds number set of simulations only (which has $\mathrm{Re}_F = 100$).

From each of the runs, we reanalyzed the data to collect instantaneous $z$-dependent profiles of the horizontally averaged flow $\bar{u}_x(z, t)$ and Reynolds stresses $\bar{R}_{ij}(z, t)$ once the system has reached a statistically stationary state. These can then be directly compared with the closure model predictions for given values of the input parameters and of the assumed eddy scale $L$.

In what follows, we use the fiducial values of the closure parameters obtained in Section 3. We then assume that $L$ is constant in the domain and find the value of $L$ that best fits a particular simulation (i.e., for a given value of $\mathrm{Ri}_F \mathrm{Pe}_F$ at $Re_F = 100$). We use a nonlinear least-squares Levenberg-Marquardt
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Values of $JPr$ measured in the DNSs of Paper I (with the solid line corresponding to the $Re_c = 9 \times 10^4$ simulations and the dashed lines corresponding to the $Re_c = 1.2 \times 10^5$ simulations) and compared with the various predicted critical thresholds for instability: the horizontal black line shows the simple model given by Equation (2) with ($JPr)_c = 0.007$; the colored lines show the predicted values of the critical threshold for stability in our closure model, namely $\sigma\epsilon_c(\sigma) C_{\nu\kappa}/C_{\nu\nu}$, for fiducial values of the parameters calibrated so far ($C_1 = 0.41$, $C_2 = 0.54$ and $C_{\nu\kappa} = 10$), for different values of $C_{\nu\nu}$. The curve with $C_{\nu\nu} = 15$ appears to be the better choice in comparison with the DNSs.</figcaption></figure>

method to iteratively select the value of $L$ that minimizes the distance between the time average of the mean flow profile $\bar{u}_x(z, t)$ obtained from the DNSs and the closure model prediction for $\bar{u}_x(z)$ once a steady state has been reached. At each iteration, we make sure to evolve the closure model equations (i.e., Equations (36)–(39) as well as Equations (42) and (44), together with the mean flow Equation (18) with $\bar{T} = 0$) in time until a steady state is reached.<fnref n="4" />

The values of $L$ thus obtained, for which the closure model is best able to predict the mean flow, are presented in Figure 7 for each simulation, as a function of $JPe_S = Ri_F Pe_F / |S|$, where $S$ is the shearing rate measured in the middle of the domain (where it is largest). We see that as expected, $L$ is commensurate with the size of the shear layer for weakly stratified flows and decreases progressively as $JPe_S$ increases. We also measured the turbulent eddy scale from the DNSs directly, using the first zero of the vertical autocorrelation function of the spanwise (y-) velocity, as in Equation (33) of Paper I. The latter is also shown in Figure 7. Remarkably, we see that as assumed in Section 3, $L$ is very well predicted by $1.5l_e$ for all available simulations. This firmly establishes that the quantity $L$ introduced in the closure model merely as a means to compute a timescale for the turbulent cascade really has a genuine physical counterpart as the vertical size of the turbulent eddies. It also implies that, moving forward, creating a good model for $L$ to be used in the closure is equivalent to creating a good model for the actual eddy scale $l_e$. This will be revisited in Section 6.

Let us now look in more detail at the comparison between the closure model prediction and the DNS results for each component of the Reynolds stress tensor. We focus on three regimes: weak ($Ri_F Pe_F = 0.01$), intermediate ($Ri_F Pe_F = 1$), and strong ($Ri_F Pe_F = 100$) stratification.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison between the best fit for $L$, the measured eddy scale $l_e$ times 1.5, and the model $L_{theor}$ (see Equation (78) in Section 5), as a function of $JPe_S$.</figcaption></figure>

## 4.1. The Nearly Unstratified Limit

Figure 8 shows in blue snapshots in time of $\bar{u}_x(z, t)$, and of the Reynolds stress components $\bar{R}_{xz}(\partial_z \bar{u}_x, t)$, $\bar{R}_{xx}(z, t)$, and $\bar{R}_{zz}(z, t)$, at selected times $t$ once the system has reached a statistically stationary state, from a DNS with $Re_F = 100$ and $Ri_F Pe_F = 0.01$ from Garaud & Kulenthirarajah (2016). Also shown by the solid black lines are the corresponding predictions from the closure model once the latter has reached a steady state, with fiducial parameters, for $L = 2.7$ (which provides the best fit for this simulation). In all that follows, we denote the steady-state limit of a quantity as predicted by the closure model with a superscript $\infty$, as in $q^\infty(z) = \lim_{t \to \infty} q(z, t)$.

As reported by Garaud & Kulenthirarajah (2016), these low stratification DNSs favor large-scale perturbations that are highly variable in time. The mean flow is nevertheless close to being sinusoidal and in phase with the forcing, which suggests that the turbulence in this limit simply acts as a turbulent diffusivity. We also see that $\bar{R}_{zz}$ is minimal in regions of low shear and maximal in regions of strong shear, showing that the perturbations are shear driven, as expected. However, the fact that these quantities are nonzero in regions where the shear is exactly null (here, at $z = \pi/2$ and $z = 3\pi/2$) also shows that there is significant turbulent overshoot from the adjacent sheared regions (also see Paper II).

The value of $L$ selected ensures that the closure model best fits the measured mean flow profile, so it is not surprising to find that the two are in very good agreement. The agreement between the model and the data for $\bar{R}_{xz}(z)$ (not shown) is similarly good, which is not surprising since $\bar{R}_{xz}(z)$ in this body-forced setup is strongly constrained independently of any other quantity. Indeed, starting from the mean streamwise equation in a steady state, we have

$$ \partial_z \bar{R}_{xz} = \frac{1}{\mathrm{Re}_F} \partial_{zz}^2 \bar{u}_x + \sin(z). \qquad (73) $$

In the limit of large Reynolds number, the diffusive term can be neglected, leaving a balance between momentum transport by Reynolds stresses and the imposed forcing. This balance can be integrated straightforwardly to yield

$$ \bar{R}_{xz} \approx -\cos(z), \qquad (74) $$

<footnote n="4">⁴ We attempted to use the nonlinear least-squares method directly on the steady-state model equations, but were not able to find converged solutions in that case, presumably due to the nonlinearity of the closure model equations.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between snapshots of the horizontally averaged horizontal flow $\bar{u}_x(z, t)$ and of the Reynolds stress components $\bar{R}_{xz}(\partial_z\bar{u}_x, t)$, $\bar{R}_{xx}(z, t)$, and $\bar{R}_{zz}(z, t)$ with steady-state predictions from the closure model $\bar{u}_x^\infty(z)$, $\bar{R}_{xx}^\infty(\partial_z\bar{u}_x^\infty)$, $\bar{R}_{zz}^\infty(z)$, and $\bar{R}_{zz}^\infty(z)$ for $\text{Re}_F = 100$, $\text{Ri}_F\text{Pe}_F = 0.01$, and $L = 2.7$, together with fiducial model parameters. The case with $C_{\text{turb}} = 0$ corresponds to the model discussed prior to Section 4.4, while the case with $C_{\text{turb}} = 0.2$ corresponds to the model discussed in Section 4.4.</figcaption></figure>

(where the integration constant can be set to 0 to ensure that $\bar{R}_{xz}$ changes sign when the shear does), so it is not surprising that this is found to be the case for both the data and the closure model. A more sensitive test of the quality of the closure model consists in plotting $\bar{R}_{xz}$ against the local shearing rate $\partial_z\bar{u}_x$. We see that in this weakly stratified system, the two are in relatively good agreement, with the closure model predictions within the limits of variability of the DNSs. However, the model exhibits a somewhat more pronounced s-shaped curve than the DNSs (for which $\bar{R}_{xz}$ is very close to being a linear function of $\partial_z\bar{u}_x$). Finally, the agreement between the predicted stationary state profiles for $\bar{R}_{xx}$, $\bar{R}_{yy}$ (not shown), and $\bar{R}_{zz}$ and the data is fairly poor. $\bar{R}_{zz}$ is significantly underestimated by the model for all $z$. The model also predicts a much stronger variation of $\bar{R}_{xx}$ with $z$ between regions of low and strong shear than what is actually observed in the DNSs, where $\bar{R}_{xx}$ is roughly constant.

## 4.2. The Stratified Limit

Figure 9 shows the same quantities as Figure 8 for a somewhat more strongly stratified case with $\text{Re}_F = 100$ and $\text{Ri}_F\text{Pe}_F = 1$. The value of $L$ for the closure model that best fits the data in this case is $L = 0.9$. The DNSs exhibit significantly

different dynamics from the more weakly stratified case. The intrinsic variability of the system in time is now much smaller, which could have been expected since increasing the stratification stabilizes the perturbations. The mean flow remains in phase with the forcing, but is now somewhat more triangular, as pointed out by Garaud & Kulenthirarajah (2016). The energy in the perturbations (as measured by $\bar{R}_{xx}$ and $\bar{R}_{zz}$) drops significantly but again not completely in shear-free regions. The closure model remains very good at predicting the shape of the mean flow $\bar{u}_x$, and is also somewhat more successful at predicting both $\bar{R}_{xx}$ and $\bar{R}_{zz}$ than in the weakly stratified case. It still globally underestimates $\bar{R}_{zz}$, however, and still underestimates both $\bar{R}_{xx}$ and $\bar{R}_{zz}$ in the regions of no shear. Finally, the predicted variation of $\bar{R}_{xz}$ with the shearing rate $\partial_z\bar{u}_x$ is also much more nonlinear than in the DNSs, with the greatest discrepancy between the two taking place in regions of weak shear.

## 4.3. The Strongly Stratified Limit

Figure 10 shows the same quantities as Figure 8 for a strongly stratified case with $\text{Re}_F = 100$ and $\text{Ri}_F\text{Pe}_F = 100$. The value of $L$ for the closure model that best fits the data in this case is $L = 0.34$. The intrinsic variability of the DNSs is now fairly minimal. As discussed by Garaud & Kulenthirarajah (2016), the mean flow is very nearly triangular, with regions of
<figure><img src="image_11.png" /><figcaption>Figure 9. As in Figure 8, but for Re<sub>F</sub> = 100, Ri<sub>F</sub>Pe<sub>F</sub> = 1, and L = 0.9.</figcaption></figure>

nearly constant shear separated by very thin layers with no shear. The energy in the perturbations (as measured through $\bar{R}_{xx}$ and $\bar{R}_{zz}$) is close (but not exactly equal) to zero where the shear is null. The closure model is now overall quite good at predicting the DNS data. The predicted mean flow has the correct triangular shape, although its amplitude is slightly overestimated by the model. The model correctly captures the overall variation of both $\bar{R}_{xx}$ and $\bar{R}_{zz}$ with z, but continues to underestimate $\bar{R}_{zz}$ somewhat everywhere. The model also predicts a nonmonotonic behavior for $\bar{R}_{xz}$ ($\partial_z \bar{u}_x$), which is not seen in the data.

### 4.4. Adding a Turbulent Diffusivity

Overall, we have found that this basic closure model performs adequately, though not as well as we had hoped for. It can reproduce the most important features of the data in the strongly stratified limit, but not so well in the weakly stratified limit. In particular, it has a tendency to underestimate the smoothness of the Reynolds stress profiles, as if it were missing a turbulent diffusion term that would act on stresses themselves. Indeed, the evolution equations for $\bar{R}_{ij}$ each contain a term in $\text{Re}^{-1} \partial_z^2 \bar{R}_{ij}$, which accounts for the microscopic diffusion of $\bar{R}_{ij}$. This term is essentially negligible for a large-enough Reynolds number, however, and therefore fails to account for the nonlocal effects we are trying to capture. For this reason, we propose to add a new term to the closure equations for $\bar{R}_{ij}$ that takes the form of a turbulent diffusivity,

with a diffusion coefficient proportional to $L \bar{R}^{1/2}$ put forward purely on dimensional grounds. With this term, the Reynolds stress evolution equation becomes

$$ 
\begin{aligned}
& \partial_t \bar{R}_{ij} + \bar{u}_k \partial_k \bar{R}_{ij} + \bar{R}_{ik} \partial_k \bar{u}_j + \bar{R}_{jk} \partial_k \bar{u}_i - \text{RiPe}(\bar{f}_i \delta_{jz} + \bar{f}_j \delta_{iz}) \\
& \qquad - \left(\frac{1}{\text{Re}} + C_{\text{turb}} \bar{R}^{1/2} L\right) \partial_{kk}^2 \bar{R}_{ij} \\
& = -\frac{C_1 \bar{R}^{1/2}}{L} \bar{R}_{ij} - \frac{C_2 \bar{R}^{1/2}}{L} \left(\bar{R}_{ij} - \frac{1}{3} \bar{R} \delta_{ij}\right) - \frac{C_\nu \bar{R}_{ij}}{\text{Re} L^2},
\end{aligned}
\quad (75) $$

where we have introduced a new model constant $C_{\text{turb}}$, which is of order unity. Because the new term is only active for nonhomogeneous shear flows, none of the results obtained regarding the local properties of the model discussed in Section 3 are affected.

We have found that, without refitting L, the value $C_{\text{turb}} = 0.2$ provides a significantly better fit to the DNS results for the stratified ($\text{Ri}_F \text{Pe}_F = 1$) and strongly stratified ($\text{Ri}_F \text{Pe}_F = 100$) runs. The model continues to underestimate $\bar{R}_{zz}$ but now better captures the presence of mixing in the regions of no shear due to turbulent overshoot. For instance, the model predictions are much closer to the profile of $\bar{R}_{xz}$ ($\partial_z \bar{u}_x$) measured in the DNSs especially in the region of no shear. The same is true for $\bar{R}_{xx}(z)$ and $\bar{R}_{zz}(z)$. On the other hand, the agreement between the model and the data is now significantly worse for the run with $\text{Ri}_F \text{Pe}_F = 0.01$, where a much better fit would be obtained with a much smaller value of
<figure><img src="image_12.png" /><figcaption>Figure 10. As in Figure 10, but for Re<sub>F</sub> = 100, Ri<sub>F</sub>Pe<sub>F</sub> = 100, and L = 0.34.</figcaption></figure>

$C_{\text{turb}} \approx 0.01$. Of course, we could improve the agreement by refitting $L$ for these low $\text{Ri}_F\text{Pe}_F$ runs, but this would invalidate our argument that $L$ is proportional to the eddy scale $l_e$, and would defeat the purpose of having a physically motivated model with reliable predictive power. Instead, we have simply demonstrated that the new model, with $C_{\text{turb}} = 0.2$, performs fairly well for $\text{Ri}_F\text{Pe}_F \ge 1$, but should not be used for lower values of $\text{Ri}_F\text{Pe}_F$.

5. Comparison of the Time-dependent Closure Model Predictions with DNSs of Stratified Periodic Shear Flows

Having constrained the model parameters and gained some experience using the proposed closure to reproduce the dynamics of very simple shear flows such as stratified plane Couette flows and Kolmogorov flows in a statistically stationary state, we now investigate its ability to model more complex problems such as the simulations presented in Paper II. As in Garaud & Kulenthirarajah (2016), Paper II studies the nonlinear evolution of a body-forced diffusive shear flow by integrating Equation (72) in a triply periodic domain. This time, the nondimensional forcing function $\hat{F}(z)$ is chosen in such a way as to generate regions of strong shear adjacent to regions of weak shear, and would drive the following laminar flow in the absence of turbulent mixing:

$$U_L(z) = \text{Re}_F \frac{\tanh\left[a \sin\left(\frac{2\pi z}{L_z}\right)\right]}{\tanh(a)}, \quad \text{where } L_z = \frac{2\pi a}{\tanh(a)}, \quad (76)$$

and where $a$ is a shape parameter. Note how the domain size $L_z$ varies with $a$ in this model. As discussed in Paper II, $U_L(z)$ tends to the Kolmogorov flow $\text{Re}_F \sin(z)$ with $L_z = 2\pi$ for $a \to 0$. As $a$ increases, the midpoint shear and the total amplitude of the flow remain the same, but increasingly large regions with little to no shear appear on either side of the shear layer. As $a \to \infty$ and $L_z \to \infty$, $U_L(z)$ tends to the hyperbolic tangent profile $\text{Re}_F \tanh(z)$ instead (see Figure 1 of Paper II). The main finding of the DNSs of Paper II was that the weakly sheared regions undergo some form of mixing even if they locally do not satisfy Zahn’s criterion $J\text{Pr} < 0.007$. The extent to which they are mixed strongly depends on the stratification —a result that prompted us to develop this closure model in the first place.

In what follows, we compare the closure model predictions against one particular DNS reported in Paper II, which has $\text{Re}_F = 100$, $\text{Ri}_F\text{Pe}_F = 10$, and $a = 2$, and which starts from initial conditions that are constructed from the laminar solution $U_L(z)$ plus small perturbations (referred to as Simulation A in Paper II). Rather than focusing on the ultimate statistically stationary state, we look at the full spatiotemporal evolution of the system. Both the DNS and the closure model simulation are initialized to be very close to the equilibrium laminar flow $U_L(z)$ plus small perturbations of similar amplitudes (note that the perturbations are added to the full 3D velocity field in the DNSs, while they are added to the horizontally averaged
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of the DNSs with the closure model. In both panels, the solid line corresponds to the DNSs, the dotted line to the closure model predictions without turbulent diffusion for the Reynolds stresses (i.e., with $C_{\text{turb}} = 0$), and the dashed line to the closure model predictions with turbulent diffusion (with $C_{\text{turb}} = 0.2$). The left panel shows the evolution of $\langle u_x^2 \rangle$ (blue), $\langle u_y^2 \rangle$ (green), and $\langle u_z^2 \rangle$ (red) from $t = 0$ to $t = 100$, while the right panel shows corresponding quantities at early times.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Comparison of the DNSs with the closure model for $C_{\text{turb}} = 0.2$ and $L = 0.7$, at times $t = 10, t = 20$, and $t = 30$. In each panel, the solid line corresponds to the DNSs, and the dashed line corresponds to the model. In the bottom-left panel, $\bar{R}_{xz}$ is shown as a function of $\partial_z \bar{u}_x$, although the axes only correspond to the data for $t = 10$. The data for $t = 20$ and $t = 30$ are shifted down and to the left by 2 and 4 units, respectively, to avoid overlap of the lines. The cross marks the shifted origin. The model captures the spatiotemporal evolution of the data within 20%–30%, except for $\bar{R}_{zz}$ (and $\bar{R}_{yy}$), where it underestimates it by a factor of about 2.</figcaption></figure>

Reynolds stress components $\bar{R}_{ij}$ in the closure model, so they cannot be exactly comparable).

We apply the closure model with the added turbulent diffusion term for the Reynolds stresses introduced in Section 4.4, and evolve the mean flow Equation (18) with

$\bar{T} = 0$, the Reynolds stress Equations (75), and the flux Equations (42) and (44) in time from the initial conditions described above. We adopt the fiducial parameters, and take $L = 0.7$ constant throughout the domain. This particular value was not fitted, but instead was selected because this was the
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Evolution of $L$ as a function of time according to the new model given by Equation (80). The left panel shows $L(z)$ at three selected times. The right panel shows $L$ as a function of $JPe_s$ at the same three times. $L$ agrees overall with $L_{theor}$ (see Equation (78)), but deviates from it in regions of very low shear, where it is controlled by turbulent diffusion from nearby regions.</figcaption></figure>

one determined to be the best fit for the corresponding Kolmogorov flow simulation with the same parameters (i.e., the same values of $Ri_F Pe_F = 10$ and $Re_F = 100$) in the previous section, which has the same laminar midpoint shear and uses the same nondimensionalization. Two cases are considered: $C_{turb} = 0$ (i.e., without turbulent diffusion of the Reynolds stresses) and $C_{turb} = 0.2$ (with it).

Figure 11 shows a comparison of volume-averaged $\langle u_x^2 \rangle$, $\langle u_y^2 \rangle$, and $\langle u_z^2 \rangle$ observed in the DNSs, with the corresponding quantities extracted from the closure model, namely

$$ \begin{aligned} \langle u_x^2 \rangle &= \frac{1}{L_z} \int_0^{L_z} (\bar{u}_x^2 + \bar{R}_{xx}) dz, \quad \langle u_y^2 \rangle = \frac{1}{L_z} \int_0^{L_z} \bar{R}_{yy} dz, \\ \langle u_z^2 \rangle &= \frac{1}{L_z} \int_0^{L_z} \bar{R}_{zz} dz. \end{aligned} \quad (77) $$

We see that the $C_{turb} = 0$ and $C_{turb} = 0.2$ cases are qualitatively similar to one another and only differ in the details (see below). Both correctly capture many (but not all) of the observed trends from the DNSs. For instance, they are fairly successful in capturing the overall evolution of the kinetic energy of both streamwise and spanwise fluid motions, at least qualitatively. However, they fail to account for a first burst of activity in the vertical motions, associated with a first decrease in the energy of streamwise motions. This feature is seen in the 3D DNSs, but not in either of the closure model runs. Paper II attributes this first burst in the DNSs to a purely 2D instability (with motion in the $x-z$ plane that is invariant in the spanwise direction), which only later becomes unstable to 3D perturbations—this is a well-known feature of the transition to turbulence in stratified shear flows (Peltier & Caulfield 2003). Since the closure model on the other hand intrinsically disallows purely 2D motion (because of the assumption of isotropization), it only properly models the 3D phase of the instability, but misses the initial 2D phase. Even so, the time at which the instability becomes fully developed and its saturation amplitude are captured relatively well by the model.

Looking at this comparison from a more quantitative point of view, we see that the closure model with $C_{turb} = 0.2$ performs somewhat better than the case with $C_{turb} = 0$ on longer timescales. In particular, the $C_{turb} = 0$ case fails to reach a

steady state and undergoes instead a series of quasi-periodic oscillations that are not present in the DNSs. Inspection of the results shows that these oscillations are due to the $C_{turb} = 0$ model being unable to pin down the edges of the laminar region, which periodically move toward and away from the midpoint of the shear layer with time. This does not happen in the $C_{turb} = 0.2$ case, presumably because the regions of low shear are never fully quiescent.

From here on we only consider the case with turbulent diffusion for the Reynolds stresses, with $C_{turb} = 0.2$, since its behavior is more regular. To see how well that model captures the vertical profiles of the mean flow and of the various Reynolds stress tensor components, we compare them to one another at times $t = 10$ (a little after the main turbulent mixing event has taken place), $t = 20$ (when the system is starting to relax to a statistically stationary state), and $t = 30$ (where it is close to being in a statistically stationary state). We see that the closure model captures the general properties of the solutions relatively well (i.e., within ~30%) except as usual for $\bar{R}_{zz}$ (and $\bar{R}_{yy}$, not shown), which is underestimated by a factor of about 2.

To summarize, we find that the closure model thus far performs adequately, reproducing the spatiotemporal evolution of both $\bar{u}_x$ and of $\bar{R}_{ij}$ typically within 20%-30% accuracy, except for $\bar{R}_{yy}$ and $\bar{R}_{zz}$, which it systematically underestimates by a factor of about 2.

## 6. How to Select $L$

Up until this point, the closure model has relied on the user to manually input the eddy scale $L$, which is not practical for astrophysical applications where one would rather automate the process. Furthermore, this scale was assumed to be constant in the domain and in time, which may not necessarily be the case if the shear varies significantly. To address both issues, we now propose a simple prognostic equation for $L$. Note that in this section we use the same nondimensionalization as in Sections 4 and 5, but the results are easy to generalize to other nondimensionalizations.

Fitting the data from Figure 7 (which shows the values of $L$ that best fit the sinusoidally forced simulations of Garaud &
Kulenthirarajah 2016), we find that L is well approximated by

$$L_{\text{theor}} = \frac{1}{l_{\text{max}}^{-1} + L_{Z}^{-1}} \quad \text{where } L_{Z} = (J\text{Pe}_{S})^{-1/2}, \qquad (78)$$

where $l_{max}$ is the size of the shear layer (which was $\pi$ for the simulations shown in Figure 7) and $J\text{Pe}_S = \text{Ri}_F\text{Pe}_F/|S|$ is the local Richardson-Péclet number based on the local shearing rate S (written in the nondimensionalization appropriate for Figure 7). This formula guarantees that $L_{\text{theor}} \to l_{max}$ when stratification is negligible and that $L \to L_Z$ when stratification is dominant. The quantity $L_Z$ should by definition be equal to 1.5 times the Zahn scale, which is equal to $l_Z = \sqrt{\frac{(J\text{Pe})_c}{J\text{Pe}_S}}$ when written nondimensionally so

$$L_Z = 1.5 \sqrt{\frac{(J\text{Pe})_c}{J\text{Pe}_S}}, \qquad (79)$$

where $(J\text{Pe})_c$ is a constant of order unity. Fitting this expression to the data from Figure 7, we find that the result coincidentally happens to take the very simple form $L_Z \approx (J\text{Pe}_S)^{-1/2}$, which implies that $(J\text{Pe})_c \approx 0.44$ (see also Paper II). Since the nondimensionalization is the same for the DNSs of Garaud & Kulenthirarajah (2016; from which Figure 7 was created) and the simulations of Paper II, we can expect that $L_Z$ and $L_{\text{theor}}$ will be given by Equation (78).

Equation (78) can be used in a number of different ways. For instance, if a typical value of the actual shearing rate of the turbulent flow is known a priori, and if the system is close to being in a statistically stationary state, then one can simply set L to be constant and equal to the corresponding $L_{\text{theor}}$ at every point in space and for all times. This is essentially what we have been doing until this point (even though the system was not always close to being statistically stationary). A more sophisticated alternative would be to let $L_{\text{theor}}$ evolve in response to the spatiotemporal variation of the shearing rate. One can also let $l_{max}$ vary with position in this expression, which would be the case for instance if the system was wall bounded (in which case $l_{max}$ would be the distance to the nearest wall). It is worth noting however that the scale $L_{\text{theor}}$ thus created is very sensitive to the value of the local shear S, which poses two problems. First, this sensitivity tends to be a destabilizing factor in any numerical scheme selected for the evolution of the closure model, and second, using $L = L_{\text{theor}}$ is somewhat unrealistic since in practice one would not expect L to vary on a length scale significantly shorter than itself.

We therefore propose as a third alternative the following evolution equation for L:

$$\partial_t L = -\frac{\bar{R}^{1/2}}{L}(L - L_{\text{theor}}) + \left(\frac{1}{\text{Re}} + C_{\text{turb}}\bar{R}^{1/2}L\right)\partial_{zz}^2 L, \quad (80)$$

where $\text{Re} = \text{Re}_F$ when applying the model in the nondimensionalization specific to this section. The first term drives L toward $L_{\text{theor}}$ on the eddy turnover time $L/\bar{R}^{1/2}$, while the second term is a turbulent diffusion term (where the $1/\text{Re}$ is there merely to regularize the diffusion coefficient in the limit where $\bar{R}^{1/2} = 0$), which smoothes out the spatial variability of L. Note that we could have introduced additional model constants in front of the first term, and we could have selected a different constant instead of keeping $C_{\text{turb}}$ in the second term. However, we have checked that these would not substantially

change the outcome, so we prefer to leave the model as simple as possible.

The complete proposed closure equations are therefore the mean flow Equation (18), the Reynolds stress Equations (75), and the flux Equations (42) and (44) together with Equation (80) for L. We evolve them in time with the same initial conditions as in the constant L case, taking in addition $L(z, 0) = L_{\text{theor}}(z)$ at $t = 0$. The results are nearly identical to those presented in Figures 11 and 12, and are not worth showing: the model performs identically well (or poorly, depending on one's perspective) in reproducing the data. As such, Equation (80) is not meant to improve the accuracy of the model, but rather, to give a simple, stable, and automatic way of predicting L. Figure 13 shows L at times $t = 10, t = 20$, and $t = 30$ as a function of z (left) and as a function of $J\text{Pe}_S$ (right). This second plot shows that L thus constructed is roughly equal to $L_{\text{theor}}$ everywhere in the domain, except in regions of very low shear (high $J\text{Pe}_S$), where it is nearly constant. This can be understood by noting that when $S \to 0$, the evolution of L is dominated by turbulent diffusion from nearby regions.

It is worth noting that in this particular simulation where the stratification is fairly large, the selected value of $l_{max}$ (which describes the intrinsic length scale associated with the shear) is irrelevant. For more weakly stratified cases, the user must specify $l_{max}$ and could select it either by considering a shear or forcing scale height (for body-forced cases), or the domain size and/or distance to the nearest boundary (for wall-bounded cases). However, since we have also shown that the model does not perform particularly well for low stratification anyway, the uncertainty in choosing $l_{max}$ in this limit is perhaps not a particularly pressing problem.

# 7. Conclusion

## 7.1. Summary

In this series of papers (Papers I and II and this work), as well as in previous works (Garaud et al. 2015; Garaud & Kulenthirarajah 2016), we have provided an in-depth investigation of turbulent mixing (of both chemical species and momentum) driven by the diffusive shear instability, excluding the effects of chemical stratification, rotation, and magnetic fields. The diffusive shear instability only exists in the low Prandtl number environment inherent to stars and gas giant planets, since it requires thermal diffusion to be strong (low Péclet number) while viscous diffusion is weak (high Reynolds number). Using DNSs, we were able to test the nonlinear instability criterion (2) and local mixing model (5) proposed by Zahn (1974, 1992). We found both to be correct in the context of numerical experiments specifically designed to reproduce the conditions under which Zahn's model should apply, namely (i) that the Péclet number be small (Garaud & Kulenthirarajah 2016) and (ii) that the turbulence properties are inherently local with a turbulent eddy scale that is much smaller than the system scale or the shear length scale (see Paper I and the works of Prat & Lignières 2013, 2014; Prat et al. 2016). On the other hand, we also showed that the model fails when this locality condition fails (Papers I and II), which happens when the eddy scale becomes large compared with the shear length scale or the system scale. This can be problematic in two cases: (a) when the stratification is weak, in which case the intrinsic turbulent eddy scale (which is theoretically given by the Zahn scale) becomes very large, or (b) in modeling mixing near and
beyond the edge of the theoretically unstable region, when the eddy scale is commensurate with the distance to that edge. In Paper II, we also showed that the nonlinear nature of the instability in the more strongly stratified case must be treated with care, as it gives rise to hysteresis—a possibility that cannot be accounted for in Zahn’s model.

In this paper, we have therefore put forward a more realistic model for mixing by diffusive shear instabilities that can (at least in principle) self-consistently address all of these problems. This model is based on the second-order turbulence closure paradigm first proposed by Ogilvie (2003) and later applied to unstratified shear flows by Garaud & Ogilvie (2005). It is expressed nondimensionally, using a unit velocity $U_0$, a unit length $L_0$, and depends on three nondimensional parameters, namely the system-scale Richardson, Reynolds, and Péclet numbers Ri, Re, and Pe given in Equations (15)–(17). We have derived the closure equations and applied the asymptotic limit of low Péclet numbers to obtain a first set of model equations given in Equations (36)–(44). By comparing the model predictions with DNSs from Paper II, we noted the need for a turbulent diffusivity term in the Reynolds stress equations and added it. When applied to plane-parallel shear flows (with the mean flow in the horizontal x-direction and shear in the vertical z-direction), one can extract a reduced set of coupled equations:

$$
\begin{align}
\partial_t \bar{u}_x + \partial_z \bar{R}_{xz} &= \frac{1}{\mathrm{Re}} \partial_{zz}^2 \bar{u}_x + (\hat{F}(z)), \nonumber \\
\partial_t \bar{R} + 2\bar{R}_{xz} \partial_z \bar{u}_x - 2\mathrm{RiPe} \bar{f}_z &= \left(\frac{1}{\mathrm{Re}} + C_{\mathrm{turb}} \bar{R}^{1/2} L\right) \partial_{zz}^2 \bar{R} \nonumber \\
&\quad - \frac{C_1 \bar{R}^{3/2}}{L} - \frac{C_\nu \bar{R}}{\mathrm{Re} L^2}, \nonumber \\
\partial_t \bar{R}_{zz} - 2\mathrm{RiPe} \bar{f}_z &= \left(\frac{1}{\mathrm{Re}} + C_{\mathrm{turb}} \bar{R}^{1/2} L\right) \partial_{zz}^2 \bar{R}_{zz} \nonumber \\
&\quad - (C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{zz} + \frac{C_2 \bar{R}^{3/2}}{3L} - \frac{C_\nu \bar{R}_{zz}}{\mathrm{Re} L^2}, \nonumber \\
\partial_t \bar{R}_{xz} + \bar{R}_{zz} \partial_z \bar{u}_x - \mathrm{RiPe} \bar{f}_x &= \left(\frac{1}{\mathrm{Re}} + C_{\mathrm{turb}} \bar{R}^{1/2} L\right) \partial_{zz}^2 \bar{R}_{xz} \nonumber \\
&\quad - (C_1 + C_2) \frac{\bar{R}^{1/2}}{L} \bar{R}_{xz} - \frac{C_\nu \bar{R}_{xz}}{\mathrm{Re} L^2}, \nonumber \\
\bar{R}_{xz} - \frac{1}{2} \partial_{zz}^2 \bar{f}_x &= -\frac{1}{2L^2} C_{\nu\kappa} \bar{f}_x, \nonumber \\
\bar{R}_{zz} - \frac{1}{2} \partial_{zz}^2 \bar{f}_z &= -\frac{1}{2L^2} C_{\nu\kappa} \bar{f}_z. \tag{81}
\end{align}
$$

for the evolution of the mean flow $\bar{u}_x$ (where the overbar is a horizontal average), the components of the Reynolds stress tensors $\bar{R}_{zz}$ and $\bar{R}_{xz}$, (twice) the turbulent kinetic energy $\bar{R} = \bar{R}_{xx} + \bar{R}_{yy} + \bar{R}_{zz}$, and the temperature fluxes $f_x = \bar{F}_x/\mathrm{Pe}$ and $f_z = \bar{F}_z/\mathrm{Pe}$. Other components of the stress tensor can be obtained if desired, but do not have to be computed otherwise. The unstratified limit ($\mathrm{RiPe} = 0$) with $C_{\mathrm{turb}} = 0$ recovers the model of Garaud & Ogilvie (2005). The model can easily be exported to a spherical coordinate system if desired, to study shear-induced mixing in stars.

The set of Equations (81) contains five closure parameters $C_1$, $C_2$, $C_\nu$, $C_{\nu\kappa}$, and $C_{\mathrm{turb}}$ and relies on the user to propose a model for the turbulent eddy scale $L$. In their local form (i.e., ignoring all spatial derivatives except $\partial_z \bar{u}_x$, which is assumed to

be constant) and assuming steady state, the closure model equations are algebraic and can be solved analytically or semi-analytically. The properties of this local model have been analyzed and recover some of the known properties of diffusive shear instabilities, including the fact that $J\mathrm{Pe}_L$ should be smaller than a constant of order unity for instabilities to proceed (where $J$ is the gradient Richardson number and $\mathrm{Pe}_L$ is the eddy-scale Péclet number; Townsend 1958; Dudis 1974; Zahn 1974). It also emphasizes the crucial role of the parameter $J\mathrm{Pr}$ (where $\mathrm{Pr}$ is the Prandtl number) in controlling the dynamics of the system, as in Zahn’s model (Zahn 1974).

Comparison of the local closure model predictions with DNSs from Paper I shows very good agreement and can be used to fit the first four model constants. We therefore propose that

$$C_1 \approx 0.41, \quad C_2 \approx 0.54, \quad C_{\nu\kappa} \approx 10 \text{ and } C_\nu \approx 15, \quad (82)$$

with errors of order 20% for $C_1$, $C_2$, and $C_\nu$, and of order 50% for $C_{\nu\kappa}$. These errors were estimated not only from comparison with the data from Paper I, but also by comparison with the constant values obtained from fitting the closure model against entirely different kinds of numerical and laboratory experiments, e.g., pipe flow and Taylor–Couette flows (Garaud & Ogilvie 2005) and Rayleigh–Bénard convection (Garaud et al. 2010).

By comparing the full closure model predictions with DNSs from Garaud & Kulenthirarajah (2016), we also found that the best fit is obtained when the eddy scale $L$ is proportional to the actual turbulent eddy scale $l_e$, measured using the vertical autocorrelation of the spanwise flow. In Paper I, we discovered that $l_e$ is directly related to the Zahn scale in the strongly stratified limit, and to the system scale in the weakly stratified limit. This then suggests that a simple theoretical model for $L$ can generally be constructed as follows:

$$L_{\mathrm{theor}} = \frac{1}{l_{\mathrm{max}}^{-1} + (1.5l_Z)^{-1}}, \quad (83)$$

where $l_{\mathrm{max}}$ is the outer scale (i.e., the overall size of the shear layer, or the shear length scale, or the distance to the wall in wall-bounded shear flows), and $l_Z = ((J\mathrm{Pe})_c|S|/\mathrm{RiPe})^{1/2}$ is the Zahn scale (where $S = d\bar{u}_x/dz$ is the local shearing rate and $(J\mathrm{Pe})_c$ is a constant of order unity), both expressed in units of $L_0$. The constant $(J\mathrm{Pe})_c$ can be fitted to the data if data are available as we have done in Section 6 for one particular set of DNSs, or taken to be equal to 1 otherwise.

As proposed, the set of Equations (81) together with Equations (82) and (83), with $C_{\mathrm{turb}} = 0$, can adequately model the dynamics of diffusive shear flows under a wide range of parameters (with the caveats discussed below) and for different kinds of setups. However, the model is not particularly stable numerically and is prone to unphysical small-amplitude oscillations associated with its inability to pin down the edges of the turbulent region adequately. Setting $C_{\mathrm{turb}}$ to a small value, e.g., $C_{\mathrm{turb}} = 0.2$, solves that problem without affecting the other results much as long as $\mathrm{RiPe}$ is of order unity or larger (for smaller $\mathrm{RiPe}$, a model with $C_{\mathrm{turb}} \equiv 0$ fares better). Similarly, letting $L$ evolve with time according to Equation (80) rather than merely setting $L = L_{\mathrm{theor}}$ also improves the stability of the numerical scheme, without noticeably changing the ultimate results. This is important for stellar evolution calculations, where stability is essential.
Finally, note that while we have so far ignored the possibility of chemical stratification and chemical mixing, this closure model can in fact easily be used to compute a turbulent diffusivity $D_{\text{turb}}$ for chemical species (assuming they are a passive tracer, i.e., assuming they do not contribute to the buoyancy much). Indeed, we have shown in Paper I (see also Prat et al. 2016) that $D_{\text{turb}} \simeq \nu_{\text{turb}}$, where $\nu_{\text{turb}}$ is the turbulent viscosity. Since we can compute $\nu_{\text{turb}}$ from our model as

$$ \nu_{\text{turb}} = - \frac{\partial \bar{R}_{xz}}{\partial S}, \qquad (84) $$

$D_{\text{turb}}$ can similarly be obtained.

## 7.2. Model Caveats and Future Work

The proposed model, as demonstrated in Sections 3 through 6, can be used to obtain qualitatively reliable predictions for the complex spatiotemporal behavior of diffusive stratified shear flows across a vast region of parameter space. For instance, it correctly captures the multiple evolutionary timescales associated with the flow, can correctly identify the edge of a mixed layer, as well as predict mixing beyond that edge. It also naturally captures the existence of multiple statistically stationary states, which can lead to hysteresis depending on initial conditions, a property of these shear flows that was first discussed in Paper II. It qualitatively captures the profile of the mean flow and of each component of the Reynolds stress tensor, but with varying degrees of quantitative accuracy.

In particular, the model is not particularly good in quantitatively reproducing the system dynamics in the limit of very low stratification. It does not correctly capture the partitioning of energy between the various components of the stress tensor (e.g., underpredicting $\bar{R}_{zz}$ and overpredicting $\bar{R}_{xx}$), and unless $C_{\text{turb}} = 0$, it does not accurately model the mean flow and momentum transport either. This may not be too much of a problem, as the very low stratification limit is likely not particularly prevalent in stars, but this limitation should be kept in mind by the user. When stratification increases, however, the model fares significantly better and is able to predict the mean flow within ~20%. It much more accurately predicts $\bar{R}_{xx}$ as well. Curiously, it continues to systematically underestimate the amplitude of $\bar{R}_{zz}$ by a factor of about 2, even though the shape is accurately predicted. Its predictive power for the shape of $\bar{R}_{xz}(S)$ (whose derivative is the turbulent viscosity) is relatively good as well.

The discrepancies between the model predictions and the DNSs must necessarily be due to some incomplete modeling of the triple correlations, pressure-strain correlations, and large-scale diffusive terms by the proposed closure. We have looked into further additions to the model to attempt to solve the problem. In particular, given that the pressure satisfies the (dimensional) Poisson equation

$$ \frac{1}{\rho_0} \nabla^2 p = - \nabla \cdot (\boldsymbol{u} \cdot \nabla \boldsymbol{u}) + \alpha g \frac{\partial T}{\partial z}, \qquad (85) $$

we see that the pressure-strain correlation terms $\overline{u_j \partial_i p}$ in Equation (24) will contain further triple-correlation terms in the unstratified limit (where $T=0$), but may also contain correlations between $u_j$ and $T$ in the stratified case, which had so far been ignored. We are presently looking into the possibility of adding new closure terms proportional to the heat

flux in the Reynolds stress equation. Preliminary results are encouraging and will be discussed in a future publication.

In general, however, our primary motivation is to create a closure that follows the same guiding principles as those outlined in the original work of Ogilvie (2003): sufficiently simple to be used in stellar evolution calculations, only containing closure terms with a simple physical explanation, and calibrated against DNSs. With this, we have shown its ability to predict the vertical profiles of the mean flow and of the stress tensor within a factor of order unity in the worst case scenario, but often much more precisely. While not perfect, this is still much better than the vast majority of mixing prescriptions used in stellar astrophysics, which have rarely or never been tested numerically.

Finally, and as discussed in Paper II, there are significant limitations to the applicability of this model that the user should bear in mind. First, it does not take into account the effect of composition in the buoyancy of the material and should therefore not be applied in regions where the density gradient due to composition is of the order of, or exceeds, the density gradient due to temperature stratification. Second, the validity of this model has only been demonstrated in the limit of low Péclet number. Whether any of it still holds for larger Péclet number flows remains to be determined. Third, we caution the reader against using it to model rotational shear indiscriminately, since the effect of rotation on the dynamics of diffusive shear instabilities has not been studied yet (this will be the subject of a future paper). To be precise, we believe that the model should almost certainly apply as is if the Rossby number based on the eddy scale $L$ is very large (i.e., if $R = \bar{R}^{1/2}/2\Omega L \gg 1$, where $\Omega$ is the mean nondimensional rotation rate of the region considered), but not if $R$ becomes of order unity or smaller. Since the shear $S$ is often dominated by rotational shear in stars, one should be particularly careful in checking whether this condition is met before applying the model. We are presently working to expand our understanding of stratified shear instabilities in all three of these directions.

L.K. and P.G. gratefully acknowledge funding by NSF AST-1517927. The simulations were run on the Hyades supercomputer at UCSC, purchased with an NSF MRI grant. The PADDI code used was generously provided by S. Stellmach. We thank D. Gagnier and J. Verhoeven for access to the result of their DNSs. We also warmly thank the IDEX initiative at Université Fédérale Toulouse Midi-Pyrénées (UFTMiP) for making this collaboration possible.

## References

Bischoff, T. 2013, in Proc. 2013 program in Geophysical Fluid Dynamics, ed. C. Cenedese, E. Chassignet, & S. Llewellyn-Smith (Woods Hole, MA: Woods Hole Oceanographic Inst.)

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

Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43

Endal, A. S., & Sofia, S. 1978, ApJ, 220, 279

Garaud, P., Gagnier, D., & Verhoeven, J. 2017, ApJ, 837, 133

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

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

Garaud, P., & Ogilvie, G. I. 2005, JFM, 530, 145

Garaud, P., Ogilvie, G. I., Miller, N., & Stellmach, S. 2010, MNRAS, 407, 2451

Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368

Knobloch, E. 1984, GApFD, 29, 105

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

Maeder, A. 1995, A&A, 299, 84

Maeder, A., & Meynet, G. 2000, ARA&A, 38, 143

Ogilvie, G. I. 2003, MNRAS, 340, 969
Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4

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

Pinsonneault, M. H., Kawaler, S. D., Sofia, S., & Demarque, P. 1989, ApJ, 338, 424

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

Zagarola, M. V., & Smits, A. J. 1998, JFM, 373, 33

Zahn, J.-P. 1974, in IAU Symp. 59, Stellar Instability and Evolution, ed.

P. Ledoux, A. Noels, & A. W. Rodgers (Dordrecht: Reidel), 185

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