Spanwise heterogeneous roughness, more specifically, spanwise-adjacent strips of relatively high and low roughness, is known to cause large-scale secondary flows in the vertical domain in neutral turbulent boundary layers. In this work, we study the response of secondary vortices to thermal stratification with focus on the stably stratified case. We systematically vary the Monin–Obukhov (MO) length scale from *L*/*h* = 0.3, which corresponds to relatively strong stratification, to ∞, which corresponds to the neutral condition. Here, *L* is the MO length and *h* is an outer length scale. The results show that stable stratification suppresses the vertical motions and reduces the vertical size and strength of the first off-wall secondary vortex. Moreover, the reduced vertical extent of the first off-wall vortex and the shear near its top together give rise to a second vortex, and the second to a third vortex. This leads to a stack of secondary vortices in a stably stratified boundary layer flow with decreasing strength when moving away from the wall.

## I. INTRODUCTION

The study of flow over heterogeneously rough surfaces is relevant in a variety of applications. A prominent example is the atmospheric boundary layer (ABL), whose characteristics can be significantly influenced by abrupt changes of terrain, e.g., due to variations in the vegetation or the presence of urban structures and bodies of water.^{1} Other examples of heterogeneous rough surfaces are turbine blades, which often undergo non-uniform surface degradation,^{2} and river beds, which are usually characterized by longitudinal bedforms and sand ribbons.^{3,4}

To provide a systematic understanding, the problem of heterogeneous roughness is widely studied under two canonical set-ups: flow over a surface whose aerodynamic properties undergo variations in the streamwise direction and one in which the variations occur in the spanwise direction. The former problem has been extensively studied in the past due to its relevance in triggering internal boundary layers.^{5–8} The latter problem, which has received an increased attention in the recent years, is associated with the formation of secondary flows as will be discussed in the following. The present research focuses on the latter problem.

It is known since the earlier works that the presence of rough stripes^{9,10} or riblets^{11} elongated in the streamwise direction can trigger secondary motions in channels and ducts. More recent experiments^{12–14} established the fact that any spanwise heterogeneity of wall conditions, e.g., existence of “high” and “low” roughness or regions with different riblet patterns, can give rise to large-scale spatial heterogeneity and secondary rolls in the mean flow. A result of this is the emergence of spanwise-alternating regions of momentum deficit and excess—termed low- and high-momentum pathways (LMP and HMP), respectively.^{14} Later on, Willingham *et al.*^{15} used large-eddy simulation (LES) to provide proof for the formation of roughness-induced secondary flows in ABL. Anderson *et al.*^{16} showed that the secondary flows observed in the previous works^{12,14,15} are of Prandtl’s second kind, i.e., caused by the gradients of Reynolds stress components as opposed to the secondary flows of the first kind, which are caused by centrifugal forces.^{17,18} Obviously, the former kind occurs exclusively in turbulent flows.

It is common to study the problem of spanwise-heterogeneous roughness under two general categories of “strip-type” and “ridge-type” roughness, to separate the “shear-inducing” and “height-increasing” effects of roughness.^{4,19–21} The former category refers to the presence of strips of rough (or rougher) surface along the streamwise direction without a considerable increase in the surface elevation. Idealized strip-type roughness studies—in which regions of “high” and “low” shear stress are imposed as boundary conditions to mimic varying roughness—are conducted by Willingham *et al.*^{15} using LES and by Chung *et al.*^{22} using a direct numerical simulation (DNS). Both studies identify large-scale rolling flows, in which the downwelling flow is located above the high-roughness strips. The same rotational sense of motion is experimentally observed in ducts^{9,10} and developing boundary layers^{20,23} lending support to the occurrence of strip-type behavior in different flow scenarios. A related configuration, studied in the context of superhydrophobic surfaces, is streamwise-elongated areas with free-slip and no-slip boundary conditions, which are addressed in a number of recent works.^{24–27} Furthermore, secondary flows are observed above walls with herringbone-patterned riblets;^{28–31} here, the spanwise varying property of the surface is not the value of friction coefficient but the directional properties of roughness.

One should note that an ideal strip-type roughness can be regarded realistic only at very large values of Reynolds number, e.g., in ABL, where terrain-related roughness is often much smaller than the boundary layer thickness, and indeed, the roughness may lie below the first wall-normal grid point in numerical simulations.^{15,32} At moderate to low Reynolds numbers, a clear-cut distinction may not be possible in reality as roughness generation inherently alters the elevation of the surface. In this sense, surface-resolving simulations of flow over heterogeneous roughness have been reported at moderate friction Reynolds numbers (*Re*_{τ} < 600) by Xie *et al.*^{33} and Stroh *et al.*,^{21} who investigated the effect of spanwise width and mean elevation of the rough strips, respectively. The latter authors identified a switch from the strip- to the ridge-type behavior as the rough strips increasingly protrude relative to the smooth strips.

Ridge-type roughness is also investigated extensively in the past using wind tunnel experiments,^{34–37} LES,^{38,39} and DNS.^{11,35,40} Here, the term “ridge” refers to streamwise-elongated areas of elevated surface height. Most previous studies indicate secondary rolls with the upwelling flow above the ridge. Notably, Yang and Anderson^{38} reported a reversal in the sense of motion at increased values of ridge spacing in their LES at very large Reynolds numbers. Later on, Anderson *et al.*^{41} developed a similarity solution for the averaged vorticity equation, which reproduces the dependence of the sense of motion on the ridge spacing. Medjnoon *et al.*^{36} also showed that in the case of adequately wide rectangular ridges, extra pairs of rolls emerge above the centerline of the ridge. Recently, Stroh *et al.*^{42} discussed the ridge-induced secondary flow as a viable heat transfer enhancement technique.

Previous studies provide a wealth of knowledge on the physics and properties of roughness-induced secondary flows. However, they entirely deal with neutral boundary layers, despite the fact that heat transfer takes place in several real-world applications leading to thermal stratification. ABL, which is the main motivation of the present work, is, for example, only neutrally stratified during relatively short transition periods after sunset or in windy conditions with a complete cloud cover. Consequently, a comprehensive judgment on the relevance and role of secondary flows in ABL calls for investigations of this phenomenon in the presence of buoyancy forces. In this sense, Bakhuis *et al.*^{43} had recently employed DNS and optical velocimetry techniques to study the turbulent Taylor Couette system (flow in the gap between two concentric rotating cylinders) with sandpaper stripes mounted on the inner cylinder. These authors examined the case of pure inner cylinder rotating, which given the analogy between the centrifugal and thermal instabilities^{44} is relevant in the understanding of unstably stratified boundary layers. They showed that the spacing of the stripes can be used to manipulate the size of turbulent Taylor vortices. Notably, they reported upwelling (outward) flow to always occur over the rough regions, which is opposite to the sense of motion in neutral flows. There is no report on the stably stratified boundary layers to the best of our knowledge.

In view of the above, the present work is devoted to understanding the response of secondary flows to stable thermal stratification. Stable boundary layer (SBL) usually occurs in the atmosphere at night and is featured by weaker turbulence compared to neutral and unstable boundary layers due to the destructive effect of negative buoyancy.^{45} SBL has been the focus of many numerical studies in the past using both LES (e.g., Refs. 46 and 47) and DNS (e.g., Refs. 48–50). However, no attention has been paid to the problem of heterogeneous roughness in this context. In the present work, we conduct LES in a periodic open-channel flow configuration with strip-type spanwise-heterogeneous roughness at very large Reynolds numbers. The studied geometry (strip size and spacing) and the Reynolds number are selected similar to those in Ref. 15 as a widely cited reference work. We aim to shed light on the characteristics of secondary flows in the presence of stable thermal stratification. To this end, we carry out simulations at different levels of stratification including the reference neutral case. The structure of the paper is as follows. In Sec. II, the details of LES as well as a description of the studied cases are presented. In Sec. III, the results of the simulations are shown. We finalize the paper in Sec. IV with the concluding remarks and an outlook to future work.

## II. COMPUTATIONAL SETUP

We use an in-house pseudo-spectral code for our LESs. The code was extensively utilized for thermally stratified boundary-layer flow calculations (see, e.g., Refs. 51–58), where we also validated the code for flows at both stratified and neutral conditions. It solves the filtered incompressible Navier–Stokes equations, under the Boussinesq approximation, and the filtered transport equation for the potential temperature in a half channel with height *h*,

where *t* is the time, $P\u0303$ is the filtered pressure, $U\u0303$ is the filtered velocity, $\Theta \u0303$ is the filtered potential temperature, *ρ* is the fluid density, *ν* is the kinematic viscosity, $D$ is the molecular diffusivity, Θ_{0} is the reference temperature, $\tau ij=UiUj\u0303\u2212U\u0303iU\u0303j$ is the subgrid-scale (SGS) stress, $qi=Ui\Theta \u0303\u2212U\u0303i\Theta \u0303$ is the SGS heat flux, *F*_{p} is an imposed pressure gradient in the streamwise direction and is kept constant in time, $u\tau =Fph$ is the total/bulk friction velocity, *S*_{Θ} = *u*_{τ}*θ*_{τ}/*h* is a source term for the potential temperature,^{57,59,60} where *θ*_{τ} is the total/bulk surface temperature scale, *i* = 1, 2, and 3 are the streamwise (*x*), spanwise (*y*), and wall-normal (*z*) directions, respectively, *δ*_{ij} is the Kronecker delta tensor, and the angle brackets represent a horizontal average. The SGS turbulent fluxes are modeled using the anisotropic minimum dissipation model (AMD).^{55,61,62} The flow is periodic in both the streamwise and the spanwise directions. The upper boundary conditions are a zero-stress zero-flux condition, $\u22023U\u03031|top=\u22023U\u03032|top=\u22023\Theta \u0303|top=0$, and zero vertical velocity, $U\u03033|top=0$.^{54,63,64} For the bottom surface, the wall shear stresses and the surface heat flux are specified according to Refs. 65–67,

Here, *τ*_{$w$,x/y} is the local wall-shear stress in the streamwise and the spanwise directions and *q*_{$w$,z} is the local wall-heat flux. $U\u0303\Delta $, $V\u0303\Delta $, and $\Theta \u0303\Delta $ are the spatially filtered LES velocity at a distance Δ*z* from the wall, Δ*z* is the height of the first off-wall grid, *z*_{o} is a pre-specified roughness length scale and models the effect of the surface roughness, $U\Vert =U\u0303\Delta 2+V\u0303\Delta 2$ is the wall-parallel velocity, *L* is the Monin–Obukhov (MO) length, and *ψ*_{M} and *ψ*_{H} are corrections that accounts for thermal stratification (here implemented in the classical form in Ref. 68). Θ_{s} is the bottom surface temperature, and it is kept constant in time. In our LESs, $u\tau =Fph=0.4$ m/s, *h* = 500 m, *ν* = 1.5 × 10^{−5} m^{2}/s, *Pr* = 0.7, *z*_{o,H} = 0.5 m, *z*_{o,L} = 0.005 m, and Θ_{s} = 293 K. It follows that the friction Reynolds number is *Re*_{τ} = 13 × 10^{6}. In this Reynolds number regime, the drag on the wall is only a function of *z*_{o}, and the flow is no longer affected by the Reynolds number, and therefore, one often calls the flow being at a nominally infinite Reynolds number (see, e.g., Refs. 38, 64, and 65,). In order to vary the strength of thermal stratification, the magnitude of *S*_{Θ} is specified for each simulation. Further details of the solver can be found in Refs. 56 and 57 and are not included here for brevity.

Figure 1 shows the computational domain. The size of the domain is *L*_{x} × *L*_{y} × *L*_{z} = 4*π* × 2*π* × 1(*h*), and the grid size is *N*_{x} × *N*_{y} × *N*_{z} = 256 × 128 × 128 following Refs. 16 and 38. Here, *h* ≡ 1 is the half channel height. The grid planes are staggered in the vertical direction. The surface roughness is heterogeneous in the spanwise direction. The roughness height is *z*_{o,H}/*h* = 10^{−3} on the two strips and is *z*_{o,L}/*h* = 10^{−5} everywhere else. The two strips are separated in the *y* direction by *s*_{y} = *π*.

Table I shows the thermal stratification condition of the LESs. We vary the bulk MO length from *L*/*h* = ∞, which corresponds to the neutral condition, to *L*/*h* = 0.3, which leads to relatively strong thermal stratification, covering the typical values encountered during a day/night cycle. Besides the MO length, the Richardson number (*Ri*) can also be used to quantify thermal stratification. For the flows considered in this work, *Ri* is not an independent flow controlling parameter. Its value is determined given the MO length *L* and other flow conditions. Generally, the Richardson number is an increasing function of *z*/*L*. Throughout the paper, all flow statistics we show are averaged in time for about 20 flow-through times after the flow reaches a statistically stationary state. Here, one flow through is *t* = *L*_{x}/*u*_{b}, and *u*_{b} is the bulk velocity.

## III. RESULTS

### A. Topology of secondary flow

The three-dimensional time-averaged flow is homogeneous in the *x*-direction and heterogeneous in the *y* (and obviously in the *z*) direction: Fig. 2 shows contours of the streamwise- and time-averaged longitudinal velocity deficit and vectors of the in-plane motions for cases N, L20, L2, and L0.3. The flow in L4 and L1 looks qualitatively similar to L2 and is not shown here for brevity. The velocity deficit *u*′′ is defined as the difference between the local mean velocity and its horizontally averaged value or

In the present notation, the overbar indicates averaging in time and *x*-direction and the angle brackets indicate averaging in the spanwise direction hereafter. For the sake of brevity, we drop the tilde sign for filtered quantities hereafter. The mean flow repeats itself from *y*/*h* = *π* to *y*/*h* = 2*π*, and we show only the flow from *y*/*h* = 0 to *y*/*h* = *π*.

Figure 2(a) compares well with the previous work (Ref. 15, case B) in terms of the sense of motion and the position of LMP and HMP: the downwelling motion and the HMP both occur above the rough strip and vice versa. This kind of behavior is generally expected when it comes to the strip-type roughness.^{9,22} In the present case, the absolute values of velocity deficit can reach up to 10% of *U*_{0} (i.e., the horizontally averaged mean velocity at *z*/*h* = 1) near the center of the HMP and LMP for the neutral flow.

The spanwise heterogeneity in $U\xaf$ is visibly a decreasing function of thermal stratification. In case L20, where the flow is weakly stratified, it is observed that the topology of the secondary motion and the alternating pattern of the HMPs and LMPs is similar to the neutral case. However, the eyeball center of the rolls is shifted toward the wall. Simultaneously, the LMP becomes weaker and its associated velocity deficit forms two distinct minima. The heterogeneity of the streamwise velocity keeps becoming less pronounced as the MO length further decreases. In case L0.3, where the flow is strongly stratified, one cannot find clear-cut HMPs and LMPs (by eyeballing the streamwise velocity contour).

The more interesting phenomenon is the stacking of the secondary vortices as the flow becomes thermally stratified. In cases N and L20, the vertical domain contains one secondary vortex, and it extends to *z*/*h* = 1. In case L2, the vertical domain contains two secondary vortices: the first extends from the wall to about *z*/*h* = 0.4, and the second extends from about *z*/*h* = 0.4 to *z*/*h* = 1. In case L0.3, if we look closely at the alternating pattern of the spanwise velocity and count, we would see three secondary vortices in the vertical domain: their strengths decrease away from the wall and the vortices near top of the domain is barely visible. We will discuss in detail the strength of vortices in Sec. III C. In their DNS of neutral boundary layer, Stroh *et al.*^{21} have recently reported the formation of additional vortices in the vertical domain when the main vortex (the one adjacent to the wall) becomes weaker, which is arguably related to the stacking phenomenon observed here.

Figures 3 and 4 further elucidate the topology of the secondary rolls by showing their height and center locations; these quantities are determined based on the minima and maxima in the profiles of wall-normal dispersive stress, which are a manifestation of the strength of the mean vertical motion. Details will be explained in Subsection III B. Figure 3 shows the height of the wall-adjacent secondary vortex (denoted as *h*_{1})—the strongest and arguably the most significant one in terms of transport from the surface—as a function of thermal stratification. In cases N (which is not shown in Fig. 3) and L20, the vertical domain contains one secondary vortex, and their height is 1(*h*). Thermal stratification limits vertical turbulent motions,^{69} and consequently, the size of the wall adjacent secondary vortex decreases. In cases L4, L2, and L1, the height of the wall adjacent vortex is about 0.5(±0.1), and in case L0.3, the height the wall adjacent vortex is 1/3. The majority of previous works on strip-type roughness report the vortex to reach the edge of the boundary layer unless when the roughness wavelength is much smaller than the boundary layer height, in which case the vortex scales down with this wavelength.^{20,22} The present results suggest that thermal stratification can act as an additional factor in limiting the vertical scale of the secondary motion and also lead to stacking of the vortices. As mentioned before, the shearing motion at the top of the wall-adjacent vortex gives rise to a second vortex, and the second to the third. The stacking of the vortices, however, is not a gradual process. From Fig. 4, we see that the center location of the top vortex remains almost constant for 1 < *L*/*h* < 4, and abruptly branches to two vortices. Here, *c*_{i} is the center location of the i-th eddy.

### B. Spanwise-averaged mean statistics

Figures 5(a) and 5(b) show the inner-scaled mean velocity and mean potential temperature profiles. Flow stratification suppresses turbulence mixing and, as a result, the mean flow (in plus units) increases as a function of flow stratification. This is still the case here. Thermal stratification affects the wall region *z* > *CL* more than the wall layer, where *C* is a *O*(0.1) constant. Hence, both the velocity and the temperature profiles collapse in the wall layer. The difference between a stratified case and the neutral case becomes large as thermal stratification gets stronger. Case L4 follows case N up to about *z*/*h* = 0.2. Case L0.3 follows case N up to only about *z*/*h* = 0.04. The effect of stratification is also parameterized by the gradient Richardson number

A large *Ri* suggests a strong effect of thermal stratification and a small *Ri* suggests a weak effect of thermal stratification. Figure 5(c) shows *Ri* as a function of *z*/*L*. *Ri* is an increasing function of *z*/*L*. This is consistent with the discussion above. From Fig. 5(c), we also observe that L0.3 behaves differently than L1–L4. While in our LES, L0.3 is fully turbulent, the large *Ri* in L0.3, particularly in the upper part of the half channel, can well lead to relaminarization of the flow at a lower Reynolds number.^{70}

In Fig. 6, we show the profiles of wall-normal momentum fluxes, that is, the spanwise-averaged Reynolds stress component $\u27e8u\u2032w\u2032\xaf\u27e9$ and the dispersive stress component ⟨*u*$\u2033$*$w$*$\u2033$⟩. Here, $w\u2033$ is defined in a similar way to Eq. (3), and $ui\u2032=Ui\u2212Ui\xaf$ are random velocity fluctuations. Note that we only show the resolved stresses and not the SGS stresses. From Fig. 6, we see that in cases N and L20, the dispersive stress is comparable to the Reynolds stress in the bulk, which highlights the role of secondary motions in momentum transport. In cases L2 and L0.3, however, the Reynolds stress dominates the dispersive stress everywhere. In fact, in cases L2 and L0.3, the dispersive stress is practically negligible for *z*/*h* > 0.5 and *z*/*h* > 0.3, respectively. These positions mark the wall-normal extent of the main (wall-adjacent) secondary vortices in each case. This is an indication that the second (and the third) secondary vortex pairs are virtually insignificant in momentum transport.

As we see in Sec. III A, stable stratification gives rise to a stack of secondary vortices, whose numbers, center locations, and strengths vary as a function of the thermal stratification. While “eyeballing” the results gives the center locations of these secondary vortices as well as their extents, it is certainly more desirable if we could resort to a more quantitative measure (that agrees with our “eyeballs”). Here, we take L2 as an example and examine a few basic statistics. Inherently, secondary flow is accompanied by a variation of the mean flow in the spanwise direction. Hence, the Reynolds stress, which measures a flow’s unsteadiness, would not be a good (direct) measure of secondary flow, and we will gain more insight if we focus on the dispersive stress. Figure 7 shows the square root of ⟨*u*$\u2033$*u*$\u2033$⟩, ⟨*v*$\u2033$*v*$\u2033$⟩, and ⟨*$w$*$\u2033$*$w$*$\u2033$⟩ in case L2. Compared to the streamwise component, the spanwise and the wall-normal components are obviously more direct measures of the secondary flow. In addition, comparing the spanwise and the wall-normal components, the wall-normal component ⟨*$w$*$\u2033$*$w$*$\u2033$⟩ appears to be a more convenient measure of the secondary flow: the locations of the two maxima measure the center locations of the two pairs of secondary vortices; the values of the two maxima measure the strength of the two secondary vortices; last, the distances between the minima, i.e., from 0 to *z*/*h* = 0.4 (the first minimum), and from *z*/*h* = 0.4 to 1, measure the wall-normal extents of the two secondary vortices. It is worth noting that defining, e.g., the size of a vortex will (always) involve some ambiguity. In the present work, we adopt ⟨$w\u2033w\u2033$⟩ as one reasonable and convenient measure, while there are certainly other options.

Figure 8 shows the profiles of ⟨$w\u2033w\u2033$⟩ for all cases. The formation of the second (in all but N and L20 cases) and the third (in case L0.3) vortex pairs in the vertical domain is evident from this figure. The location and extent of these vortices are already discussed in Sec. III A. It is also worth noting that the ⟨$w\u2033w\u2033$⟩ profiles collapse in the wall layer, similar to the mean velocity profile. We also observe that the peak values of the dispersive stress are smaller for cases with stronger stratification and also as distance from the wall increases in a stack of vertices. This peak value can be considered as a measure of the strength of secondary flow. In Subsection III C, we will study in detail the effect of stable stratification on the strength of the secondary flow.

### C. Strength of secondary vortices

At the neural condition, the domain contains one pair of secondary vortices between two roughness strips. If the distance between the two roughness strips are not too close,^{20,22} the secondary vortices extend to the top of the computational domain, the centers of the secondary vortices usually locate close to *z*/*h* = 0.5, and integrating the absolute mean streamwise vorticity (or compensated circulation^{36,38}) over the spanwise and wall-normal coordinates gives a fairly good measure of the vortices’ strength. Hence, measuring secondary vortices is less of a problem at the neural condition. For stably stratified flow, because a secondary vortex does not occupy the whole vertical height of the half channel, if one was to measure their strength by integrating the vorticity field, one would have to manually cut the domain in order to isolate each vortex, which could be cumbersome. As a result, we adopt the peak values of the square root of dispersive stress (see Fig. 8) as the measure of strength of the secondary motion in this paper. It might be indeed the most appropriate choice as the dispersive stress is directly linked to the vertical motion and the resulting momentum, heat, or mass transport.

Figure 9 shows the strengths of the secondary vortices in all stratified cases normalized by that of the neutral case. For flows that have more than one secondary vortex in its vertical domain, the first off-wall vortex is, by far, the most prominent. Thermal stratification weakens the first off-wall vortex, and its strength is a monotonically increasing function of the MO length. According to the PI theorem, the strength of the first off-wall vortex is a function of thermal stratification only, or equivalently a function of *L*/*h* only, for a specified domain and specified roughness at sufficiently high Reynolds numbers. We get 0.42(*L*/*h*)^{0.16}, i.e., a weak power law dependence, by fitting our limited data. The MO length *L* varies by about *O*(100) during typical day/night cycles. According to our fitting, this translates to a factor of 2 in terms of the secondary flow’s strength. The variation in the strength of the second off-wall vortex is also depicted as a function of stable stratification in Fig. 9. Here, we observe that the strength of the second off-wall vortex is not necessarily monotonic unlike the first off-wall vortex. Among the studied cases, *L*/*h* = 2 generates the strongest second vortex.

### D. Structure of turbulence

From the previous discussions, we see that thermal stratification suppresses the secondary motion above the spanwise heterogeneous roughness, which also suppresses the three dimensionality of the resulting turbulent Reynolds stress. To see that more clearly, we show the invariant map for cases N and L0.3. Here, $\eta =(aijaji/6)1/2$ and $\xi =(aijajkaki/6)1/3$ are the invariants of the anisotropic part of the normalized Reynolds stress tensor, *a*_{ij}, as defined in Ref. 71. Besides *η* and *ξ*, *II* = −*a*_{ij}*a*_{ji}/2 and *III* = *a*_{ij}*a*_{jk}*a*_{ki}/3 are also the invariants of *a*_{ij}.^{72} We see from Fig. 10 that case N shows much three dimensionality in the Reynolds stress, whereas L0.3 is almost identical to that of a plane channel at a neutral condition.^{71} It can be seen that the structure of turbulence is particularly different in the bulk region toward the edge of the boundary layer where it approaches the one-component turbulence for the stable case.

Figure 10 is based on the spanwise-averaged Reynolds stress components, i.e., $\u27e8ui\u2032uj\u2032\xaf\u27e9$. In order to shed light on the local structure of turbulence, we show the invariant maps based on the local Reynolds stresses (i.e., $ui\u2032uj\u2032\xaf$) at fixed lateral [(a), (b), (d), and (e)] and wall-normal [(c) and (f)] locations in Fig. 11. Comparing the vertical profiles for the neutral case [(a) and (b)] with those of the stable case [(d) and (e)] confirms that the thermal stratification leads to more one-dimensional turbulence far from the wall. On the other hand, when comparing the two lateral positions for the neutral case, it appears that the one above the smooth strip tends more toward the axisymmetric two-component turbulence at the edge of boundary layer. From Figs. 11(d) and 11(e), the structure of the turbulence looks to be more similar at different *y*-locations. To further clarify this point, we plot Lumley’s map for all the points at the *z*-location of the center of the main secondary vortex, which is displayed in Figs. 11(c) and 11(f). The distribution of points in these plots clearly shows a more uniform structure in the stable flow (f) compared to the neutral one (c).

## IV. CONCLUSIONS

We report the LES results of a neutral and stably stratified turbulent flow above alternating strips of relatively low and high surface roughness elongated in the streamwise direction. At the neutral condition, as shown by previous researchers, the spanwise heterogeneity of roughness gives rise to secondary vortices and the resultant HMPs and LMPs. The secondary vortices extend to the edge of the boundary layer—meaning that there is only one vortex in the vertical domain—and contribute to the vertical momentum transport across the boundary layer.

Our results generally show that stable stratification tends to reduce both the strength and the size of the secondary vortices but does not alter their sense of motion. At a weakly stratified condition (*L*/*h* = 20), the topology of in-plane flow is relatively similar to that of the neutral condition. As stratification intensifies, it reduces the vertical size of the wall-adjacent vortex and gives rise to one (at *L*/*h* = 1, 2, 4) or two (at *L*/*h* = 0.3) additional secondary vortices in the vertical domain forming a “stack” of vortices. Notably, in all these cases, the vertical sizes of all vortices in the stack are about the same. For *L*/*h* = 1, 2, and 4, the vertical size of the two secondary vortices is about *h*/2, and for *L*/*h* = 0.3, the size is about *h*/3. Moreover, in the two-vortex regime, the wall-normal locations of the vortex centers are nearly insensitive to stratification.

Despite the formation of the vortex stack due to stratification, the first off-wall vortex remains by far the strongest vortex in the stack; it’s strength (measured by the square root of the peak value of wall-normal dispersive stress) is typically an order of magnitude larger than that of the second vortex. When there is a third vortex, its strength is smaller than the second one, but their ratio is not as larger as that of the first and the second. It seems physically justified for the wall-adjacent vortex to be the strongest as it is directly driven by the wall heterogeneity, while the other vortices form as a result of the shear on top of the first one. Importantly, it is only the first off-wall vortex that causes a significant vertical momentum transport. This implies that the stacking phenomenon can lead to a step change in the wall-normal extent, by which secondary motions can contribute to the transport phenomena. That is up to half and one-third of the boundary layer thickness in two- and three-vortex regimes, respectively.

The strength of the first off-wall secondary vortex is a decreasing function of thermal stratification (or an increasing function of *L*/*h*). For the roughness considered here, i.e., *s*_{y} = *πh*, *z*_{o,H}/*h* = 10^{−3}, and *z*_{o,L}/*h* = 10^{−5}, it is a function of only *L*/*h*. The data suggest that the strength follows ∼(*L*/*h*)^{0.16}, i.e., the strength halves if *L* increases by *O*(10^{2}) (a typical variation during a day/night cycle). This, along with the changes in the topology discussed above, signifies meaningful modifications in the characteristics of the secondary flow in ABL. Nevertheless, the rather weak dependence means that the secondary flow remains relevant up to rather high levels of stable stratification.

Finally, one should note that the present work is a first attempt to shed light on the effect of thermal stratification on roughness-induced secondary flows. This is arguably an essential first step to understand the role of secondary flows in many real-world problems, most notably in ABL. Motivated by this particular application, we adopted LES at infinitely large Reynolds number and focused on a strip-type configuration. This, obviously, does not rule out the need for further investigations (e.g., DNS at low Reynolds numbers and wind tunnel measurements at intermediate Reynolds numbers) to get a better picture of the problem.

## ACKNOWLEDGMENTS

X.Y. thanks ONR (Grant No. N000142012315) for the financial support.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.