**2. LES governing equations**

Along with laminar-turbulent transition, turbulence parameterization constitutes the most critical part of flow modeling [68, 80]. In high-Reynolds-number turbulent flows computational limitations impose the choice of a grid size <sup>∆</sup>*grid* substantially larger than the smallest scale of motion. In LES, the separation of scales between resolved and unresolved scales is achieved by filtering (with a spatial filter of characteristic size <sup>∆</sup> <sup>∆</sup>*grid*) the equations describing the transport of mass, momentum and scalar quantities. In particular, the filtered equations (using the Boussinesq approximation) governing continuity, conservation of momentum and scalar transport are

$$\frac{\partial \tilde{u}\_{l}}{\partial \mathbf{x}\_{l}} = \mathbf{0} \,, \tag{1}$$

$$\frac{\partial \widetilde{u}\_{i}}{\partial t} + \frac{\partial \widetilde{u}\_{i} \widetilde{u}\_{j}}{\partial \mathbf{x}\_{j}} = -\frac{\partial \widetilde{p}}{\partial \mathbf{x}\_{i}} - \frac{\partial \tau\_{\overline{i}j}}{\partial \mathbf{x}\_{j}} + \nu \frac{\partial^{2} \widetilde{u}\_{i}}{\partial \mathbf{x}\_{j} \partial \mathbf{x}\_{j}} + \widetilde{f}\_{i} \ . \tag{2}$$

$$\frac{\partial \widetilde{\theta}}{\partial t} + \widetilde{u}\_{\dot{l}} \frac{\partial \widetilde{\theta}}{\partial \mathbf{x}\_{\dot{l}}} = -\frac{\partial q\_{\dot{l}}}{\partial \mathbf{x}\_{\dot{l}}} + \kappa \frac{\partial^2 \widetilde{\theta}}{\partial \mathbf{x}\_{\dot{l}} \partial \mathbf{x}\_{\dot{l}}} \, \, \, \tag{3}$$

where (*<sup>u</sup>*1, *<sup>u</sup>*2, *<sup>u</sup>*<sup>3</sup>)=(*<sup>u</sup>*, *<sup>v</sup>*, *<sup>w</sup>*) are the components of the resolved velocity field, *<sup>θ</sup>* is the resolved scalar, *p* is the effective pressure, *ν* is the kinematic viscosity, *κ* is the scalar diffusivity, and *f <sup>i</sup>* is a forcing term. In the stable/unstable ABL flow, the buoyancy force and the Coriolis force would be included as *f <sup>i</sup>* = *δi*3*g θ* −�*<sup>θ</sup>* �*H <sup>θ</sup>*<sup>0</sup> <sup>+</sup> *fcεij*3*<sup>u</sup><sup>j</sup>*, where *<sup>θ</sup>* represents the resolved potential temperature, *<sup>θ</sup>*<sup>0</sup> is the reference temperature, �·�*<sup>H</sup>* denotes a horizontal average, *g* is the gravitational acceleration, *fc* is the Coriolis parameter, *δij* is the Kronecker delta, and *εijk* is the alternating unit tensor. In homogeneous rotating turbulence, the Coriolis force would be included as *f <sup>i</sup>* = −2*εij*3Ω*ju<sup>k</sup>*, and without loss of generality, we would chose −→<sup>Ω</sup> = (0, 0, <sup>Ω</sup>). The effects of the sub-grid scales on the evolution of *<sup>u</sup><sup>i</sup>* and *<sup>θ</sup>* appears in the SGS stress *τij* and the SGS flux *qi*, respectively. They are defined as

$$
\pi\_{i\underline{j}} = \widetilde{u\_i}\widetilde{u}\_{\underline{j}} - \widetilde{u}\_{\bar{i}}\widetilde{u}\_{\underline{j}\prime} \qquad \text{and} \qquad q\_{\underline{i}} = \widetilde{u\_i}\widetilde{\theta} - \widetilde{u}\_{\bar{i}}\widetilde{\theta}.\tag{4}
$$

SGS models are needed to parameterize *τij* and *qi* as a function of the resolved velocity and scalar fields.
