Turbulence Structure and Mixing in Strongly Stable Boundary-Layer Flows over Thermally Heterogeneous Surfaces

Direct numerical simulations (DNS) at bulk Reynolds number Re =104\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^4$$\end{document} and bulk Richardson number Ri =0.25\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=0.25$$\end{document} of plane Couette flow are performed with the results used to analyze the structure and mixing intensity in strongly stable boundary-layer flows. The Couette flow set-up is used as a proxy for a real-world stable boundary layer flow with surface thermal heterogeneity. Along the upper and lower walls, the temperature is either homogeneous or varies sinusoidally, but the horizontal-mean surface temperature is the same in all cases. Over homogeneous surfaces, the strong stratification always quenches turbulence resulting in linear velocity and temperature profiles. However, over a heterogeneous surface turbulence survives. Molecular diffusion and turbulence contribute to down-gradient momentum transfer. The total (diffusive plus turbulent) heat flux is directed downward, but its turbulent contribution is positive, i.e., up the mean temperature gradient. Analysis of covariances of velocity and temperature, their skewness, and the flow structure suggests that counter-gradient heat transport is due to quasi-organized cell-like vortical motions generated by surface thermal heterogeneity. These motions transfer heat upwards similar to their counterparts in highly convective boundary layers. Thus, the flow over heterogeneous surface features local convective instabilities and upward eddy heat transport, although the overall stratification remains stable with downward mean heat transfer. The DNS results are compared to the results from large-eddy simulations of weakly stable boundary layers (Mironov and Sullivan in J Atmos Sci 73:449–464, 2016). The DNS findings corroborate the key role of temperature variance in setting the structure and transport properties of stably stratified flow over heterogeneous surfaces, and the importance of third-order transport of the temperature variance.


Introduction
Modelling the stably stratified planetary boundary layer (SBL) is of considerable importance for numerical weather prediction, climate studies, and related applications. Model errors associated with the SBL, a regime characterized by relatively weak turbulence and low mixing intensity, are often substantial (Holtslag et al. 2013), and it is still largely unclear how the trouble can be cured. Many SBL features are poorly understood (see, e.g., Mahrt 2014, for discussion). One particularly challenging issue is the strongly stable boundary layer over a thermally heterogeneous surface. In the present study, we use direct numerical simulation (DNS) to gain insight into the effect of surface thermal heterogeneity on the structure and transport properties in a strongly stable boundary layer. Mironov and Sullivan (2016, hereafter MS16) used large-eddy simulation (LES) to examine turbulence in the atmospheric SBL over thermally homogeneous and thermally heterogeneous surfaces. The LES data were used to perform a comparative analysis of secondmoment budgets and mixing intensity in homogeneous and heterogeneous SBLs. A physically plausible explanation of the enhanced mixing in the heterogeneous SBL was found, and possible ways to parameterize the heterogeneity effects in atmospheric models were discussed. It should be emphasized that the results of MS16 are pertinent to weakly-to-moderately stable boundary layers characterized by continuous and rather vigorous turbulence. The structure and mixing intensity in strongly stable boundary layers characterized by weak and often intermittent turbulence are still largely unknown and need to be investigated. The present study attempts to make a step forward in this direction.
Among the many questions to explore, the following outstanding questions should be addressed: (i) If turbulence dies out over a homogeneous surface due to strong stability, does it survive over a heterogeneous surface? (ii) If turbulence survives, how anisotropic is it and does it generate appreciable vertical fluxes of momentum and scalars? (iii) What are the particular features of the second-moment budgets in a strongly stable regime, and what is the role of pressure-scrambling and third-order transport in maintaining the variance and flux budgets? The present contribution addresses issues (i) and (ii). Analysis of the second-moment budgets is left for future work.
The present study is not aimed at simulating a specific real-world geophysical (e.g., atmospheric or oceanic) flow. Instead, we focus on physical processes at work in strongly stratified boundary-layer flows, namely on the effect of surface thermal heterogeneity on the structure and transport properties of turbulence. To this end, we use an idealized plane Couette flow set-up (a physical analog of our numerical configuration is discussed in the next section). The flow is driven by a fixed velocity of the upper surface, while the lower surface is at rest. The stable density stratification is imposed by a fixed temperature difference between the upper and lower boundaries. The temperature at the horizontal upper and lower surfaces is either homogeneous or varies sinusoidally in the streamwise direction, while the horizontal-mean temperature is the same in all cases. The DNS data are used to compute vertical profiles of mean fields, second-order and (some) third-order statistical moments of turbulence, and to analyze the structure and mixing intensity over thermally homogeneous and thermally heterogeneous surfaces.
Clearly, a real-world atmospheric SBL and an idealized Couette flow differ in many respects (e.g., the surface roughness can play an important role). We tend to think, however, that results from the present analysis with respect to the maintenance of turbulence in heterogeneous strongly stable flows are relevant to the real-world flows, at least qualitatively. The relevance of our DNS findings for atmospheric SBL flows is further discussed in Sect. 5.
A few words are in order to make clear what is referred to as "turbulence" in the context of the present study. By applying averaging over horizontal planes to compute turbulence statistics (see Sect. 3), we basically assume that turbulence incorporates all fluctuations within the DNS model domain. That is, turbulence encompasses both nearly isotropic and nearly homogeneous small-scale fluctuations and quasi-organized coherent structures that are strongly non-isotropic and non-homogeneous (a discussion of coherent structures encountered in our Couette-flow configuration is given in Sect. 4). Different definitions are used by different research communities, however. Quite often, quasi-organized coherent motions would not be classified as turbulence in the atmospheric science literature. Those motions are referred to as "sub-meso-scale motions" to discriminate them from the classical Kolmogorov-type nearly isotropic and nearly homogeneous turbulence.
From the standpoint of a large-scale or meso-scale atmospheric model, for which the entire DNS domain is a single model grid box, the definition adopted in the present study would mean that turbulence comprises all sub-grid scale (SGS) fluctuations. Statistical moments of SGS fluctuating velocity and scalar fields (fluxes and variances) should be described by a turbulence parameterization scheme. Note, however, that atmospheric models with relatively coarse resolution do not apply a single parameterization scheme that describes the effects of various types of SGS fluctuations in a unified framework. A number of parameterization schemes are applied instead, where each scheme serves to describe the effects of a particular process separately, e.g., the effects of internal gravity waves, cumulus convection (shallow and deep), and turbulence. This artificial separation of processes in numerical models of the atmosphere causes many conceptual and practical problems (see, e.g., Arakawa 2004;Mironov 2009, for discussion). It is desirable to achieve a more coherent description of the various types of SGS fluctuations within a unified parameterization framework. The present-day tendency with high-resolution models is to switch off gravity-wave and cumulusconvection parameterization schemes (at least deep convection scheme), or to make those schemes less active, assuming that at high resolution the SGS motions are more isotropic and homogeneous. Then, the relative importance of turbulence parameterization schemes considerably increases. It should be emphasized, however, that even with the finest resolution attainable today (and likely for many years to come) for numerical weather prediction and climate models the SGS fluctuating velocity and scalar fields remain anisotropic and heterogeneous, and this is especially true for stably stratified boundary layers. It is this type of boundary layers that is the focus of the present study.
The present paper is closely related to the legacy of Sergej Zilitinkevich, who had keen interest in the stably stratified boundary layers throughout his scientific career. Sergej's contributions to the SBL studies include the now classical SBL depth scale (Zilitinkevich 1972), resistance and heat-transfer laws for stably stratified boundary layers (e.g., Zilitinkevich 1989a, b), and turbulence closure models for stably stratified flows (e.g., Zilitinkevich et al. 2013), to mention a few. Sergej showed particular interest in the problem of maintenance of turbulence in strongly stable boundary-layer flows and the role of coherent quasi-organized motions in setting the SBL mean and turbulent structure. In this regard (and in relation to the present paper), the work of Glazunov et al. (2019) should be mentioned, where DNS of plane Couette flow is used to get insight into the SBL structure. Now we give an outline of the paper. The flow configuration and governing equations are described in Sect. 2. Details of the simulations performed are given in Sect. 3. Simulation results are discussed in Sect. 4. Conclusions are presented in Sect. 5.

Flow Configuration and Governing Equations
We consider a three-dimensional, stably stratified, plane Couette flow. The fluid depth is H , the lower boundary is at rest, and the upper boundary moves with a constant velocity U u . Stable buoyancy stratification is maintained by a temperature difference Δθ = θ u − θ l between the upper and lower boundaries. The physical characteristics of the fluid are the kinematic molecular viscosity ν, the molecular temperature conductivity κ, and the buoyancy parameter β i = −g i α T , where g i is the acceleration due to gravity and α T is the thermal expansion coefficient. Both ν and κ are taken to be constant in our simulations. The Boussinesq approximation is used, and the simplest linear equation of state is utilized, where ρ is the fluid density, ρ r is the constant reference density, θ is the potential temperature (for the sake of brevity, it will also be referred to as simply "temperature"), and θ r is the constant reference potential temperature. Both α T and β i are constant.
The governing equations given below contain three dimensionless parameters. These are the bulk Reynolds number, the Prandtl number, and the Richardson number, where g 3 is the magnitude of the gravity vector; in the co-ordinate system used here g i = (0, 0, −g 3 ).
A real-world flow configuration closely analogous to our numerical configuration can hardly be found in geophysical flows. A conceivable physical analog of our numerical set-up is a laboratory tank (a channel-like apparatus of infinite band type) filled with fresh water at room temperature. Using H = 0.25 m, U u = 0.04 m s −1 , and ν = 10 −6 m 2 s −1 (the value for fresh water at 20 • C), we obtain Re = 10 4 . Using a quadratic fresh-water equation of state (Mironov et al. 2010 is an empirical coefficient and θ r = 277.13 K is the temperature of maximum density of fresh water, we find that a temperature difference between the upper and lower horizontal boundaries of about 1 K is needed to obtain Ri = 0.5. To arrive at this estimate, we recast the Richardson number in terms of the buoyancy difference between the upper and lower horizontal boundaries, Ri = H Δb/U 2 u , where b = g 3 (ρ r − ρ) /ρ r , g 3 = 9.81 m s −1 , θ l = 292.65 K (19.5 • C), and θ u = 293.65 K (20.5 • C). Note that the Prandtl number is about 7 for fresh water (at 20 • C), which is considerably higher than the value of Pr = 1 that we adopt for our numerical fluid.
The flow variables are made dimensionless with the length scale H , velocity scale U u , time scale H /U u , and temperature scale Δθ , and dimensionless temperatureθ = (θ − θ l ) / (θ u − θ l ) is introduced. The governing equations written in dimensionless form are (in order to simplify the notation, we omit hats over dimensionless variables): Here, t is time, x i are the right-hand Cartesian co-ordinates, u i are the velocity components, and p is the kinematic pressure (deviation of pressure from the hydrostatically balanced reference pressure divided by the constant reference density). The Einstein summation convention for repeated indices is adopted. The origin of the co-ordinate system is at the lower boundary, the x 3 axis is aligned with the gravity vector and is positive upward, and the x 1 axis is in the direction of U u . Note that, since the vertical-velocity equation in our simulations is solved for the fluctuation of u 3 about its horizontal mean, θ r can be set equal to θ l so that the dimensionless reference temperature (θ r − θ l ) / (θ u − θ l ) drops out from the buoyancy term on the right-hand side (r.h.s.) of Eq. 4. Periodic boundary conditions for u i and θ are applied in both x 1 and x 2 horizontal directions. At the horizontal boundaries, the following Dirichlet boundary conditions are used: and: where L 1 is the domain size in the x 1 direction, δθ is the (dimensionless) amplitude of the temperature variations at the upper and lower surfaces, and n is the number of cold and warm stripes (the number of surface temperature waves). In the homogeneous case, δθ = 0. In the heterogeneous cases, δθ > 0 but the horizontal-mean surface temperature is the same as in the homogeneous case. Notice the lower and upper boundary conditions (7) and (8) preserve the flow symmetry of a Couette flow configuration. Along the upper boundary the temperature waves propagate in the x 1 direction with boundary speed U u which generates a thermally heterogeneous surface matching its counterpart along the lower wall. The Couette set-up although idealized has desirable features for the present application: the average boundary temperatures are constant in time, and with no mean pressure gradient forcing the average total momentum and temperature fluxes are constant with height. As a result, we are able to acquire stationary statistics by long temporal integration even under strongly stable intermittent flow conditions. Couette flow with constant momentum and temperature fluxes is an approximate analog to the assumed constant flux layer in a high Reynolds number boundary layer flow (e.g., Wyngaard 2010, p. 215).

Simulations Performed
The DNS code and its parallelization used in the present study are described in detail in Sullivan et al. (2000), Sullivan and McWilliams (2002), and Sullivan and Patton (2011). A description of the code is not repeated here; readers are referred to the above papers. One simulation with homogeneous lower and upper surfaces (HOM) and three simulations with heterogeneous surfaces (HET) are performed. In all simulated cases, the fixed values of Pr = 1, Re = 10 4 , and Ri = 0.25 are used. The number of grid points is (512, 512, 256) in the (streamwise x 1 , spanwise x 2 , vertical x 3 ) directions, respectively. The domain size in the vertical direction is 1, and the domain size in both horizontal directions is 8. The number of cold and warm stripes (the number of surface temperature waves) in the heterogeneous cases is 4. Governing parameters of the simulations performed are summarized in Table 1.
In Table 1, T t is the total length of the simulation in dimensionless time units, T s is the length of the sampling period (at the end of the run), Re τ = u * Re is the wall Reynolds number based on the surface friction velocity u * , Ri τ = u −2 * Ri = Re −2 τ Re 2 Ri is the Richardson number based on the surface friction velocity u * (e.g., Gandía-Barberá et al. 2021), Q * is the surface temperature flux, and 1/L is a bulk stability measure where L = −u 3 * /κRiQ * is the Obukhov length (Obukhov 1954), and κ = 0.4 is the von Kármán constant. Recall that all variables are made dimensionless with the scales H , U u , H /U u , and Δθ for length, velocity, time, and temperature, respectively.
The homogeneous simulation starts with a fully developed, stationary, neutral Couette flow. The stable buoyancy stratification is established by gradually (linearly in time) increasing the Richardson number over 100 dimensionless time units from Ri = 0 to Ri = 0.25. The simulations are then continued until turbulence dies out, and a laminar Couette flow regime is achieved. The value of Ri = 0.25 proves to be sufficient to fully quench turbulence in the homogeneous case. Based on the stability measures Ri τ and 1/L, all of the simulations in Table 1 are in the very stable regime. As a comparison, Nieuwstadt (2005) finds that turbulence is not sustained when 1/L > 0.5 in channel flow: we modified the definition of L in Nieuwstadt (2005) to match the present definition.
The heterogeneous flows start with Ri = 0, δθ = 0, and a linear velocity profile. In order to assist initial turbulence spin-up, velocity and temperature fluctuations taken from the neutral turbulent Couette flow are added in the lower 1/4 and the upper 1/4 of the computational domain. The Richardson number is increased (linearly in time) from Ri = 0 to Ri = 0.25 over 10 time units, while the temperature difference δθ is increased from zero to its value given in Table 1 over 100 time units. The simulations over heterogeneous surfaces are carried forward over many time units until a quasi-stationary flow regime is reached; the flow is considered stationarity when the average momentum and scalar flux profiles are constant with height. The number of time samples differs between the cases, but the sampling period in the heterogeneous cases covers more than 150 time units (see Table 1). The DNS data are averaged over horizontal planes, and the resulting profiles are then averaged over several thousand time steps. These horizontal and time mean quantities are treated as approximations to the ensemble-mean quantities. In what follows, an overbar denotes a horizontal-mean quantity, a prime denotes a fluctuation about a horizontal mean, and the angle brackets denote quantities averaged over time.
Note that, apart from the results from the homogeneous simulation HOM and three heterogeneous simulations HET025, HET050, and HET075 (which are the major focus of the present study), results from the fully turbulent neutral and convective Couette-flow simulations are also shown for comparison in Figs. 3 and 8. The neutral run (used to initialize the stably stratified homogeneous run HOM) and the convective run (used to initialize the neutral run) are similar to HOM, but the Richardson numbers are Ri = 0 and Ri = − 0.2, respectively, and the length of the sampling period is T s = 100.

Mean Fields
Vertical profiles of the streamwise mean-velocity component U = u 1 and mean temperature Θ = θ are shown in Fig. 1. The spanwise mean-velocity component is negligibly small in all simulations and is not shown. In the homogeneous case, the profiles of both U and Θ are linear, corresponding to the laminar solution of the plane Couette problem. Although we provide turbulence a good chance to survive by starting the homogeneous simulation with a vigorously turbulent neutral flow and gradually increase the Richardson number, the buoyancy stratification at Ri = 0.25 is strong enough to fully extinguish turbulence.
As seen from Fig. 1, the mean velocity is only slightly affected by the surface thermal heterogeneity, the mean flow gradient in the middle of the domain is only slightly larger than the laminar Couette value, ∂U /∂ x 3 > 1. The situation is noticeably different for mean temperature. As the amplitude δθ of the surface temperature variations increases, the flow becomes increasingly mixed with respect to Θ away from the walls. The associated increase of the temperature gradient near the horizontal boundaries is due to the fixed temperature boundary conditions. Note that the effect of temperature heterogeneity on mean temperature is appreciable because of the relatively large values of δθ in simulations HET050 and HET075. In simulation HET025, δθ is comparatively smaller, the effect of surface heterogeneity is weak, and the temperature profile is nearly linear. Vertical profiles of the (shear, buoyancy) frequencies (S 2 , N 2 ) and gradient Richardson number Ri g , shown in Fig. 2, further expose the impact of surface heterogeneity on the mean flow state, especially in the wall region. Here, Inspection of Fig. 2 reveals a surprisingly complex dependence on δθ, not apparent in Fig. 1, especially for the shear frequency. Notice the variation of S 2 at a fixed vertical location near the boundaries below x 3 < 0.2 is not monotonic with increasing δθ, e.g., at x 3 = 0.1 S 2 = (1.02, 0.25, 0.51) with increasing amplitude of the heterogeneity. At the same height, N 2 = (0.25, 0.36, 0.36). As a result, the vertical profile of the gradient Richardson number varies in a complex way and simulation HET050 features the largest most stable Ri g ∼ 1.4 at locations x 3 = (0.1, 0.9) near the walls. The turbulent fluctuations generated near the boundary do impact the mean flow gradients in the middle of the domain, Ri g smoothly decreases as the amplitude of the heterogeneity increases. Figure 3 shows profiles of the streamwise mean-velocity component in wall units The flow in the close vicinity of the boundaries is very well resolved in our simulations. The first grid point above the surface is x + 3 ≈ 0.2. A larger grid spacing may be used in neutral Couette flows, but the resolution cannot be compro- mised in our simulations because of the need to resolve the large temperature gradients in the near-wall thermal boundary layers. It is these temperature gradients that drive turbulence in the heterogeneous simulations. Test runs with lower resolution in the vertical direction yield results (not shown) that are quantitatively and even qualitatively different from the results of our high-resolution simulations and are not trustworthy. Note the domain size in the horizontal directions cannot be compromised either because of the need to simulate largescale elongated structures characteristic of plane Couette flows (e.g., Komminaho et al. 1996;Papavassiliou and Hanratty 1997;Pirozzoli et al. 2014;Avsarkisov et al. 2014;Jiménez 2018;Lee and Moser 2018).
Black solid and dashed lines in Fig. 3 show, respectively, the streamwise velocity from the neutral Couette-flow simulation and the logarithmic velocity profile: where B 0 = 5.0 is a dimensionless constant. The neutral velocity profile closely follows classic log-layer scaling over the height range 20 ≤ x + 3 ≤ 200. This is not the case for the flows over thermally heterogeneous surfaces, where the velocity profiles reveal no log-layer scaling. The green dot-dashed curve in Fig. 3 shows the Monin-Obukhov log-linear velocity profile (Obukhov 1954;Monin and Obukhov 1954;Monin and Yaglom 1971): where C u = 5 is a dimensionless constant, and L + = L Re τ is the Obukhov length in wall units. The Monin-Obukhov velocity profile in Fig. 3 uses (u * , Q * ) from simulation HET050 and external parameters (Re, Ri). Inspection of the results in Fig. 3 shows the velocity profile in case HET050 (green solid curve) does not show similarity to (11). This is not surprising however, because the Monin-Obukhov surface-layer flux-profile relationships are only applicable to turbulent layers over homogeneous surfaces. For thermally heterogeneous surfaces, different flux-profile relationships are required.

Second-Order Moments
Vertical profiles of turbulence kinetic energy TKE = 1 2 u 2 i are shown in Fig. 4. TKE increases with increasing amplitude of the surface temperature variations. The turbulence level is very low in case HET025, as δθ = 0.25 proves to be too small to make the flow vigorously turbulent.
Velocity variances uu = u 2 1 (streamwise), vv = u 2 2 (spanwise), and ww = u 2 streamwise velocity variance. In case HET025, the dominate contribution to the TKE is uu, with only small contributions from vv and ww. Figure 6 shows the velocity variances in wall units (i.e., uu + = uu/u 2 * , vv + = vv/u 2 * , and ww + = ww/u 2 * as function of x + 3 = x 3 Re τ ). The streamwise velocity variances do not collapse on the same curve, but the results from simulations HET050 and HET075 (green and red curves in Fig. 6) get somewhat closer to each other if the surface friction velocity is used for normalization. A similar statement can be made about the vertical velocity variances (although with a less confidence), but not about the spanwise velocity variances.
In strongly stable flows, Mahrt (2014) finds an abundance of structures with intermittent turbulence mixed with wavelike motions and anisotropic two-dimensional modes. To quantify the turbulence anisotropy in our flows, we use the departure-from-isotropy tensor (see, e.g., Pope 2000): recall that all components of b i j are zero in isotropic turbulence. In the two-component limit, where velocity fluctuations in one direction (e.g., vertical) are suppressed, the corresponding diagonal component of b i j is equal to −1/3. As seen from Fig. 7, b 33 is negative throughout the flow in simulations HET050 and HET075, indicating that the vertical-velocity fluctuations are strongly damped by buoyancy and the turbulence is mainly confined in a horizontal plane. In simulation HET025, b 33 is slightly positive in small regions near the boundaries, but overall turbulence in HET025 is very weak and all three velocity variances and hence the TKE are negligibly small, see Figs. 4 and 5. An informative diagnostic used to characterize the turbulence anisotropy are the second II = b i j b i j and third III = b i j b jk b ik invariants of the departure-from-isotropy tensor (e.g., Lumley and Newman 1977). These invariants are also often defined as: II = − 1 2 b i j b i j and III = 1 3 b i j b jk b ik (e.g., Lumley 1978;Choi and Lumley 2001). Alternatively, one can use the quantities ξ = 1 6 b i j b i j 1/2 and η = 1 6 b i j b jk b ik 1/3 . An anisotropy map in the ξ −η plane is known as the Lumley triangle (Pope 2000), although one side is not straight. Figure 8 shows Blue circles indicate that turbulence in the heterogeneous case HET025 is either in a 1C or 2C state. The two-component state in HET025 is of little interest as the regions of the flow with 2C turbulence have negligibly small TKE. The turbulence in HET050 and HET075 is in a 2C state near the walls but smoothly transitions into an axisymmetric state similar to the neutral case above the boundary. Notice, in all three heterogeneous cases near the mid-plane x 3 = 1/2 3D turbulence is nearly quenched and the residual motions are 1C u 1 sloshing motions.
As expected, the temperature variance θθ = θ 2 increases with increasing amplitude of the surface temperature variations δθ, see Fig. 9. It is significant that the large values of θθ are confined to the immediate vicinity of the lower and upper boundaries, while the temperature variance is small in the bulk of the flow interior 0.1 ≤ x 3 ≤ 0.9.
Vertical profiles of the streamwise momentum-flux component (denoted wu) and vertical temperature flux (denoted wθ) are shown in Figs. 10 and 11, respectively. The spanwise momentum-flux component (not shown) is negligibly small in all simulations. In the steady state, the total fluxes, i.e., the sum of turbulent and molecular fluxes, are depth constant. That is, In the laminar Couette flow, the contributions to wu and wθ are solely due to molecular diffusion with total fluxes equal to Re −1 and (PrRe) −1 , respectively. As seen from Fig. 10, the magnitude of the total streamwise momentum-flux component increases with the increasing amplitude of the surface temperature variations. In the simulation HET025, the turbulent momentum flux is very small and the total flux is virtually the same as in the laminar flow. In the simulations HET050 and HET075, both the turbulent and The total vertical temperature flux, see Fig. 11, decreases with the increasing amplitude of the surface temperature variations. This is counter to the total streamwise momentum-flux component. The total temperature flux in simulation HET025 is very close to that in laminar flow. However, the turbulent temperature flux u 3 θ is not entirely negligible. A remarkable feature of the heterogeneous simulations is the sign of the vertical turbulent temperature flux. Although the flow is stably stratified in the mean and the total (molecular plus turbulent) vertical temperature flux is negative, the turbulent heat flux proves to be positive. Turbulent motions generated by the surface thermal heterogeneity transfer heat up the gradient of the mean temperature. It somewhat resembles convective boundary-layer flows, where quasiorganized cell-like structures induce counter-gradient heat transport. Figure 12 presents vertical profiles of vertical-velocity and temperature skewness:

Vertical-Velocity and Temperature Skewness and Flow Structure
As seen from the plots, in all three heterogeneous simulations (S w , S θ ) > 0 near the lower boundary, and by symmetry about the mid-plane x 3 = 1/2, (S w , S θ ) < 0 near the upper boundary. S w > 0 indicates that positive (upward) vertical velocity has a lower fractional area coverage (more localized) than negative (downward) vertical velocity. Likewise, S θ > 0 indicates a stronger localization of positive temperature fluctuations (about the horizontalmean temperature) as compared to negative temperature fluctuations. The statistics from our stably stratified Couette flow with surface thermal heterogeneity are broadly similar to a dry convective boundary layer (CBL) driven by surface potentialtemperature (buoyancy) flux (e.g., Lenschow et al. 2012) or Rayleigh-Bénard convection between plates of fixed temperature (e.g., Moeng and Rotunno 1990). Both S w and S θ are positive in the major part of the CBL, where highly localized positive potential-temperature anomalies are collocated with positive vertical velocity, forming familiar convective plume structures which account for most of the upward potential-temperature flux.
Flow visualization of the (u 3 , θ) fields in our heterogeneous simulations reveals strongly localized quasi-coherent flow structures characterized by large positive fluctuations in vertical velocity and temperature. As a result, these structures generate positive turbulent temperature flux in the lower part of the flow (Fig. 11), i.e., precisely where both S w and S θ are positive (Fig. 12). In the upper part of the flow, large negative vertical velocity with S w < 0 is collocated with large negative potential-temperature fluctuation S θ < 0, leading to a positive turbulent temperature flux.
The instantaneous snapshots in Figs. 13 and 14 reveal spatially intermittent turbulent flow structures in the wall region. Figure 13 shows fluctuations of vertical velocity, u 3 , and of potential temperature, θ , in an x 1 − x 2 plane just above the lower surface for simulation HET050; for reference, the streamwise variation of the surface temperature θ s f c = θ(x 3 = 0) given by (8) is also shown in the figures. The most vigorous positive and negative vertical velocity fluctuations are concentrated at the same streamwise locations x 1 ≈ (0.9, 2.9) with quiescent nearly laminar motions at intermediate locations x 1 ≈ (2, 4). Thus, the vigorous fluctuations in u 3 occur at a streamwise location roughly 70 degrees forward of the peak surface temperature but behind the location of maximum − ∂θ s f c /∂ x 1 at x 1 = (1, 3). Visualization also shows positive fluctuations in streamwise velocity u 1 positioned farther downstream than the peak in θ . Apparently horizontal advection in the wall region transports high amplitude turbulence forward of the location of the peak surface temperature variation, similar to weakly stratified heterogeneous flows, e.g., Stoll and Porté-Agel (2009) and MS16, and the temperature field responds more rapidly to the surface boundary conditions than streamwise velocity. It is important to mention that the spatially intermittent flow structures in our very stable heterogeneous simulations are noticeably different from the tilted temperature fronts found in weakly stable boundary layers and stratified Couette flows with continuous turbulence (e.g., García-Villalba et al. 2011;Chung and Matheou 2012;Sullivan et al. 2016;Glazunov et al. 2019).
Coherent patterns of strong positive u 3 nearly collocated with the patterns of strong positive θ are readily identified, and their spatial correlation results in positive (upward) turbulent temperature flux over a considerable part of the flow as shown in Fig. 14. Although positive values of u 3 θ cover a smaller fractional area (more localized) than negative values of u 3 θ , their amplitude is large. The peak positive values of u 3 θ are roughly four times larger in magnitude compared to their peak negative counterparts; the fluctuations u 3 θ are an order of magnitude larger than their horizontal average u 3 θ and the horizontal-mean turbulent temperature flux is positive (upward). Importantly, the flow remains stably stratified in a global (horizontal mean) sense, and the total, i.e., turbulent plus molecular, temperature flux remains negative (downward). The flow visualization supports the explanation of positive turbulent temperature flux suggested by the vertical-velocity and temperature skewness. The visualization in panel b) of Fig. 14 shows vertical turbulent temperature flux for simulation HET075. Although the pattern of u 3 θ differs in detail between simulations HET050 and HET075, e.g., in terms of distance between the "arrow-shaped" regions of strong positive u 3 θ and an increased amplitude in HET075, the overall picture of temperature flux remains the same between the two simulations.

Flux of Temperature Variance
Using LES, MS16 performed a comparative analysis of the second-moment budgets and mixing intensity in weakly (moderately) stable boundary layers over thermally homogeneous and thermally heterogeneous surfaces. Among other things, the study revealed a key role of the temperature variance in turbulent mixing in a horizontally heterogeneous SBL and the importance of the third-order turbulent transport term in maintaining the temperaturevariance budget. Due to the surface heterogeneity, the third-order moment, i.e., the vertical flux of temperature variance (denoted by wθθ), is nonzero at the surface. As a result, the turbulent transport term (divergence of the temperature-variance flux) not only redistributes the temperature variance in the vertical, but is a net gain.
The following expression is used to estimate the vertical flux of temperature variance on the basis of LES data: where a tilde denotes a resolved-scale (filtered) quantity, and the superscript "s" denotes a sub-grid (sub-filter) scale fluctuation. As in the previous sections, an overbar denotes a horizontal-mean quantity, a prime denotes a fluctuation about a horizontal mean, and angle brackets denote quantities averaged over time. The first two terms on the r.h.s. of Eq. 16 are zero at the surface becauseũ 3 = 0. The last term cannot be estimated from LES but is presumably small (see MS16 and Machulskaya and Mironov 2018, for discussion). The third term is zero in the homogeneous SBL becauseθ = 0 at the surface. In the heterogeneous SBL, surface temperature variations modify local stability conditions and thus modulate the surface temperature flux. The surface temperatureθ and the surface temperature flux u s 3 θ s prove to be positively correlated, leading to a positive flux of temperature variance at the surface.
Within the DNS framework, the expression for the vertical flux of temperature variance reads: It is instructive to establish the correspondence between the LES and DNS estimates of the temperature-variance flux. The first two terms on the r.h.s. of Eq. 16 describe the transport of resolved (first term) and sub-grid (second term) temperature variance by the resolved vertical velocityũ 3 . The sum of these terms is in one-to-one correspondence with the first term on the r.h.s. of Eq. 17 that describes the transport of temperature variance by the vertical velocity (there are no sub-grid scale quantities in DNS, hence there is only "total" velocity and "total" temperature variance).
There is no one-to-one matching of the third and the fourth terms on the r.h.s. of Eq. 16 with the second term on the r.h.s. of Eq. 17. It can be readily seen, however, that under certain Fig. 15 (Left) third-order vertical velocity-temperature covariance (vertical flux of the temperature variance) from simulations HET025 (blue), HET050 (green), and HET075 (red). Solid curves show total (turbulent plus molecular) covariance, and dot-dashed and dashed curves show contributions due to turbulence and due to molecular diffusion, respectively. (right) The same as in the left panel but for the lower part of the domain for simulations HET050 (green) and HET075 (red) conditions the LES and DNS terms are closely analogous. If the resolution is sufficiently high, most energy-containing scales of motion are simulated explicitly by LES with nearly isotropic sub-grid scale turbulence. Then, the sub-grid scale temperature flux is of diffusive and down-gradient character, that is, u s i θ s = −κ s ∂θ/∂x i , whereκ s is the sub-grid scale temperature diffusivity. The third term on the r.h.s. of Eq. 16 becomes: Applying a diffusive down-gradient approximation to the fourth term on the r.h.s. of Eq. 16, we obtain: If the temperature diffusivityκ s does not change considerably in space (κ s is small) the second and the third terms on the r.h.s. of Eq. 18 and the second term on the r.h.s. of Eq. 19 can be neglected. If changes ofκ s in time can also be neglected (κ s − κ s is small), then Eqs. 16, 18 and 19 yield the following expression for the LES-based vertical temperature-variance flux: Here, Pr s and Re s are the Prandtl number and the Reynolds number, respectively, defined in terms of (constant) sub-grid scale viscosity and (constant) sub-grid scale temperature diffusivity. Thus, the estimates of the vertical temperature-variance flux based on DNS, Eq. 17, and on LES, Eq. 20, coincide up to definitions of the temperature variance and of the Prandtl and Reynolds numbers.
Vertical profiles of wθθ from our heterogeneous DNS runs are shown in Fig. 15. Both the turbulent and molecular contributions to wθθ given by the first and second terms on the r.h.s. of Eq. 17, respectively, are positive close to the lower boundary; by symmetry, these contributions are negative close to the upper boundary. Importantly, the molecular temperature-variance flux is nonzero at the surface. Hence, the third-order transport term not only redistributes the temperature variance in the vertical, but is a net gain. The situation is broadly similar to that in weakly (moderately) stable flows, where the surface thermal heterogeneity produces wθθ that serves to increase the temperature variance in the boundary layer.

Conclusions
Direct numerical simulations at bulk Reynolds number Re = 10 4 and bulk Richardson number Ri = 0.25 are performed to analyze the structure and mixing intensity in strongly stable boundary-layer flows over thermally homogeneous and heterogeneous surfaces. An idealized plane Couette flow set-up is used as a proxy for real-world flows. The flow is driven by a fixed velocity at the upper surface, while the lower surface is at rest. The temperature at the horizontal upper and lower surfaces is either homogeneous or varies sinusoidally in the streamwise direction, while the horizontal-mean temperature is the same in the homogeneous and heterogeneous cases.
The stratification is strong enough to quench turbulence over homogeneous surfaces, resulting in linear velocity and temperature profiles. However, turbulence survives over heterogeneous surfaces. Both the molecular diffusion and the turbulence contribute to the downward, i.e., the down-gradient, transfer of horizontal momentum. The total (diffusive plus turbulent) heat flux is directed downward. However, the turbulent contribution to the heat flux appears to be positive, i.e., up the gradient of the mean temperature. An analysis of the second-order velocity and temperature covariances and of the vertical-velocity and temperature skewness suggests that the counter-gradient heat transport is due to quasi-organized cell-like vortex motions generated by the surface thermal heterogeneity. These motions act to transfer heat upwards similar to quasi-organized cell-like structures that transfer heat upwards in convective boundary layers. This finding is corroborated by visualization of the velocity and temperature fields. Thus, the flow over heterogeneous surface features local convective instabilities and upward eddy heat transport, although the overall stratification remains stable and the heat is transported downward in the mean.
Due to the surface thermal heterogeneity, the vertical flux of temperature variance is nonzero at the surface. As a result, the transport term (divergence of the temperature-variance flux) in the temperature-variance budget not only redistributes the temperature variance in the vertical, but is a net gain. The same effect is encountered in weakly (moderately) stable flows over thermally heterogeneous surfaces, where the third-order transport terms serves to increase the temperature variance in the boundary layer.
A few words are in order concerning the relevance of our DNS findings to the real-world atmospheric flows. For instance, the temperature contrast between the warm and cold stripes in our simulations may look unrealistically high at first glance. The following simple example demonstrates, however, that the governing parameters of our idealized flow configurations may be quite similar to the governing parameters of the real-world atmospheric SBL flows. Take the following values that are fairly typical of an atmospheric stably stratified boundary layer during calm clear nights: the SBL depth is 70 m, the wind speed at the SBL top is 4 m s −1 , and the temperature difference across the SBL (from top to bottom) is 2 K. With the gravitational acceleration of 9.81 m s −2 and the reference temperature of 280 K, we obtain Ri = 0.31, which is close to Ri = 0.25 in our simulations. With a 2 K temperature difference across the SBL, a 4 K difference between warm and cold stripes is needed to obtain the flow configuration very similar to that of our simulation HET050. Such a surface temperature difference is clearly possible; furthermore, even larger temperature differences at real-world surfaces are quite likely. In order to make a more definitive quantitative statement, data from targeted field measurements are required. Rough estimates may be obtained using crowdsourced data.
A subject of our future work is a comprehensive analysis of the second-moment budgets, where the emphasis is on the turbulence anisotropy and the role of pressure-scrambling effects in maintaining the budgets. It should also be examined whether major conclusions from the present analysis remain valid at higher values of the Reynolds number and for other flow configurations, e.g., for pressure-gradient driven channel flows over homogeneous and heterogeneous surfaces. Further open issues are the effects of the surface roughness, of the surface heterogeneity orientation (e.g., temperature-wave crests normal vs. parallel to the mean flow), and of the size and form of the surface patches on the flow structure and mixing intensity. Finally, efforts should be made to use the DNS findings to improve parameterizations of strongly stable boundary layers in large-scale atmospheric models.