Direct numerical simulation of stratified Ekman layers over a periodic rough surface

Abstract Ekman layers over a rough surface are studied using direct numerical simulation. The roughness takes the form of periodic two-dimensional bumps whose non-dimensional amplitude is fixed at a small value ($h^+=15$) and whose mean slope is gentle. The neutral Ekman layer is subjected to a stabilizing cooling flux for approximately one inertial period ($2{\rm \pi} /f$) to impose the stratification. The Ekman boundary layer is in a transitionally rough regime and, without stratification, the effect of roughness is found to be mild in contrast to the stratified case. Roughness, whose effect increases with the slope of the bumps, changes the boundary layer qualitatively from the very stable (Mahrt, Theor. Comput. Fluid Dyn., vol. 11, issue 3–4, 1998, pp. 263–279) regime, which has a strong thermal inversion and a pronounced low-level jet, in the flat case to the stable regime, which has a weaker thermal inversion and stronger surface-layer turbulence, in the rough cases. The flat case exhibits initial collapse of turbulence which eventually recovers, albeit with inertial oscillations in turbulent kinetic energy. The roughness elements interrupt the initial collapse of turbulence. In the quasi-steady state, the thickness of the turbulent stress profiles and of the near-surface region with subcritical gradient Richardson number increase in the rough cases. Analysis of the turbulent kinetic energy (TKE) budget shows that, in the surface layer, roughness counteracts the stability-induced reduction of TKE production. The flow component, coherent with the surface undulations, is extracted by a triple decomposition, and leads to a dispersive component of near-surface turbulent fluxes. The significance of the dispersive component increases in the stratified cases.

and, under typical sunny conditions, the air's potential temperature decreases with height. At night, however, this effect is reversed and surface cooling leads to the formation of a stably stratified inversion layer which can suppress turbulent motions. Similar suppression of ABL turbulence is seen in winter-time conditions and polar regions. The conditions under which surface cooling can completely or intermittently relaminarize the flow remain the subject of current research.
The meteorology of the stable boundary layer includes states of continuous turbulence and intermittent turbulence (e.g. Businger 1973). Mahrt (1998) categorized two prototypical states: the weakly stable boundary layer and the very stable boundary layer. The very stable boundary layer is characterized by weak winds and clear skies that lead to strong net radiative cooling at the surface. The weakly stable boundary layer, which has weaker net radiative cooling, is described by Monin-Obukhov similarity theory (Monin 1970), in which turbulence, although reduced in strength, occurs continuously in time. On the contrary, the very stable boundary layer is characterized by global intermittency where turbulence is reduced for periods which are long compared with the time scale of individual eddies (Mahrt 1989). Internal gravity waves are also present in the very stable boundary layer.
The turbulent boundary layer over a rough surface has been reviewed by Raupach, Antonia & Rajagopalan (1991) and Jiménez (2004), among others. One measure of roughness elements is their height as quantified by a roughness Reynolds number, h + = h/δ ν , the ratio of the height of the elements to the viscous scale, δ ν . With increasing h + , the boundary layer changes from a smooth to a transitionally rough to a fully rough boundary layer. The fully rough boundary layer has enhanced momentum transport and drag. For a given h + , the shape and spatial distribution of the roughness elements are important as reviewed by Flack & Schultz (2010) among others. In the present work we consider a periodic array of sinusoidal bumps with moderate slope (figure 1). The analogous problem of unstratified flow over a wavy bottom has received attention in experiments (Zilker, Cook & Hanratty 1977;Gong, Taylor & Dornbrack 1996) and, more recently, in numerical studies using DNS or wall-resolved large-eddy simulation (LES) (e.g. De Angelis, Lombardi & Banerjee 1997;Sullivan, McWilliams & Moeng 2000;Napoli, Armenio & De Marchis 2008;Yang & Shen 2010). Turbulent flow over rough surfaces has a component in the time-averaged field which is coherent with the surface structure and gives rise to a dispersive component of the turbulent fluxes. The coherent velocity can be extracted using a triple decomposition, and the so-called dispersive component is a significant contributor to turbulent fluxes in the near-surface layer (Finnigan 2000;Sullivan et al. 2000;Poggi, Katul & Anderson 2004;Li & Bou-Zeid 2019).
There is not much systematic study of the competing effects of stable stratification and roughness in canonical problems. Ohya, Neff & Meroney (1997) and Ohya (2001) performed laboratory experiments of stratified boundary layers that develop in a wind tunnel over smooth and rough surfaces, respectively. The bottom wall in the rough case had a two-dimensional roughness imposed by a chain of oval rings and had a colder temperature than the bulk flow with θ varying between 27.4 K and 44.1 K. Vertical profiles showed reduction of turbulence levels with increasing stability in both rough and smooth cases. Sullivan & McWilliams (2002) conducted DNS of a turbulent Couette flow over waves (including the stationary-wave case) under moderate stable and unstable stratification. Their results show a decrease (increase) of turbulence levels under stable (unstable) stratification.
Turbulence in the boundary layer can collapse in the presence of sufficiently strong stability. Banta et al. (2007) found long periods of suppressed turbulence during nights with a strongly stable ABL. The near-surface layer appears to decouple from the Here L x , L y and L z are the domain size in the streamwise, spanwise and vertical directions, respectively. The wavelength of the harmonic function that generates the bump is λ = L x /N b , where N b is the number of bumps. The half-length of the bump is l = λ/4. upper regions of the ABL during these episodes of suppressed turbulence. The bulk Richardson number (Ri b , defined later by (2.11)) is a critical parameter; large values of Ri b are indicative of strong stability. During periods with Ri b > 2 in the cooperative atmosphere-surface exchange study in 1999 (CASES-99) observational campaign, the coupling between the near-surface layer and the outer layer was found to be weak (Poulos et al. 2002). Two layers of turbulent kinetic energy (TKE), one above and the other below a local minimum of TKE, were identified by Banta et al. (2007) and Cuxart & Jiménez (2007). Such a two-layer configuration of TKE was also found in a recent DNS of the stratified Ekman boundary layer by Gohari & Sarkar (2018).
The so-called Ekman boundary layer (EBL) is a simplified example of the ABL, whereby a boundary layer in a rotating reference frame develops under unidirectional horizontal flow in geostrophic balance (Coriolis acceleration is equal and opposite to the pressure gradient, both being orthogonal to the flow). The stable EBL is a canonical problem for studying the stabilization of the ABL. Direct numerical simulation studies of the stable EBL have imposed buoyancy with a constant temperature boundary condition (Ansorge & Mellado 2014;Deusebio et al. 2014;Shah & Bou-Zeid 2014) or with a constant cooling flux (Gohari & Sarkar 2017; these studies were performed with a smooth-bottom boundary. The boundary-layer response to stability in these studies spanned various regimes of turbulence, depending on the relative strength of the stability: initial decrease of turbulence and even collapse of turbulence to a laminar state; recovery to continuous turbulence; recovery to global intermittency where turbulent near-surface patches co-exist with laminar flow; and, finally, episodes of complete turbulence collapse followed by recovery to spatially intermittent turbulence. Stratified turbulent channel flow, a canonical problem to investigate buoyancy effects in wall-bounded flows, has received much attention (Garg et al. 2000;Armenio & Sarkar 2002;Nieuwstadt 2005;Flores & Riley 2011;García-Villalba & del Álamo 2011;He & Basu 2015). Armenio & Sarkar (2002) show initial collapse of turbulence in a stratified channel flow followed by resurgence of turbulence. They find an outer layer with suppressed turbulence and wavy motion where the gradient Richardson number (Ri g defined by (2.8)) is larger than 0.2, and an inner layer with active turbulence where Ri g decreases from 0.2 to a small value at the wall. García-Villalba & del Álamo (2011), in their large-domain simulations, find global intermittency with laminar patches interspersed within turbulence when the stratification is strong. The surface cooling flux can be used to define the Obukhov length (L, defined later by (2.6)), and its value (L + , defined later by (2.7)) relative to the viscous scale is an important measure of the strength of buoyancy. Flores & Riley (2011) propose a criterion for turbulence collapse based on L + decreasing to 100 during the initial transient. This occurred in their stratified channel with initial L + = 683. Coleman, Ferziger & Spalart (1990), Ansorge & Mellado (2014), Shah & Bou-Zeid (2014), Deusebio et al. (2014) studied the stratified Ekman layer using DNS with a constant temperature imposed at the wall. In this problem the surface buoyancy flux decreases with time although Ri b is constant. Intriguing spatial intermittency was observed by Ansorge & Mellado (2014) in their DNS study of a stably stratified Ekman layer. Similar to Nieuwstadt (2005) and Flores & Riley (2011), initial turbulent collapse followed by recovery was observed when the surface temperature of a neutral Ekman layer was suddenly dropped to impose the prescribed value of stability, Ri b . The spatial characteristics of intermittent turbulence were then analysed in detail during this transient process. In neutrally stratified Ekman layers, Deusebio et al. (2014) found large-scale roll structures with one dominant frequency that matched the convective frequency of the low-level jet (LLJ). They also found that these counter-rotating streamwise vortices influence the near-wall structures by pushing or lifting fluid close to the wall. Shah & Bou-Zeid (2014), from analysis of the TKE budget, show that the reduction of turbulence levels in the stratified EBL is primarily due to the inhibition of shear production rather than the buoyant TKE destruction. Interestingly, this feature of buoyancy-induced reduction of turbulence production, which is key to the suppression of turbulence by density stratification is a generic feature of stratified shear flows and has been found by us in uniform shear flow (Jacobitz, Sarkar & VanAtta 1997) and later in the shear layer (Brucker & Sarkar 2007) and the stratified wake (Brucker & Sarkar 2010). The decrease in turbulence production with increasing stratification is an indirect effect of buoyancy that decreases the correlation coefficient between streamwise and vertical velocity fluctuations. Gohari & Sarkar (2017) performed DNS of the Ekman layer with a constant buoyancy flux, and found differences of this constant-flux stability case with respect to previous cases with constant temperature stability: the low-level jet is stronger, there are recurring episodes of collapse and rebirth of turbulence during the transient, and the TKE profile has local peaks at two vertical locations.
Recently, Gohari & Sarkar (2018) conducted DNS of a smooth-surface EBL that is subject to a finite-time (approximately, one inertial period) cooling flux. They found that initial L + cri = Lu * /ν 700 provides a cooling flux that is sufficiently strong to cause the initial collapse of turbulence independent of Reynolds number, Re * , where L is the Obukhov length scale and u * is the friction velocity. The turbulence collapse criterion for stratified Ekman flow is considered based on the normalized Obukhov length scale (Obukhov 1971), L + , and has similar values as inferred from the observational study of Banta et al. (2007) and the DNS of Flores & Riley (2011). The final state, for a fixed L + and a fixed cooling flux, was found to depend on Re * , because an increase in initial Re * (under the constraint of fixed L + ) is equivalent to an increase in Ri b . In particular, an EBL with a final stability of Ri b ≥ 2 relaminarized.
None of the Ekman layer DNS have considered surface roughness, which is often a feature of the ABL. This motivates the present research that addresses how the destabilizing effect of surface roughness competes with the stabilizing effect of surface cooling in the EBL. The simulation results are related to meteorological characteristics that are distinctive features of the stable ABL: low-level jets (Smedman, Tjernström & Högström 1993;Cuxart et al. 2000;Banta et al. 2002), collapse of surface-layer turbulence (Banta et al. 2007), local turbulence peaks at locations above the surface layer (Mahrt 1985;Smedman et al. 1993;Banta et al. 2007) and unsteadiness of turbulence statistics (Banta 2008;Pichugina et al. 2008;Sun et al. 2012). These characteristics not only change the atmospheric dispersion of tracers and pollutants, as has been well documented in the past, but also, as discussed in recent work, strongly influence acoustic propagation in the nocturnal boundary layer (Talmadge et al. 2008) with this influence being modified by hilly terrain (Damiens, Millet & Lott 2018).
The work is organized as follows. Section 2 describes the governing equations and the problem set-up. Section 3 describes the results of the simulations through the evolution of mean and turbulence statistics. Section 4 describes the flow structures by isosurface visualizations and elucidates the role of coherent structures by a triple decomposition. Finally, § 5 contains the discussion and conclusions.

Formulation
In the study of the stratified EBL, several important physical parameters arise: the geostrophic wind U ∞ , the Coriolis frequency f , the turbulent Ekman layer thickness δ * and the friction Reynolds number Re * . The governing equations for the conservation of momentum under the Boussinesq approximation and potential temperature in a rotating reference frame are as follows: Here t is time, x j is the spatial coordinate in the j direction, u j is the velocity component in that direction, p is the pressure deviation from the mean pressure imposed by geostrophic and hydrostatic balance, δ i3 is the Kronecker delta, ij3 is the alternating unit tensor, ν is the molecular viscosity, β is the thermal expansion coefficient for air, g is the gravitational acceleration, f is the Coriolis parameter, α is the thermal diffusivity and θ is the deviation of potential temperature from its constant reference value. It can be shown that Reynolds number and Prandtl number, defined as are the only non-dimensional parameters controlling the dynamics of a neutral Ekman flow at steady state, when the statistics have been adequately decorrelated from their initial condition (Spalart 1988). Although the laminar Ekman layer depth, D, is not a proper length scale describing the momentum transport in a turbulent Ekman flow, Re D provides a universal comparison point among different Ekman flow studies. For a turbulent Ekman flow, it is proper to use the turbulent Ekman layer thickness, δ N = u * /f , where u * is the friction velocity which is defined below and subscript N denotes neutral conditions. The friction velocity and the friction Reynolds number (Re * ) are computed as Hereafter, · denotes the Reynolds average, which is computed as a horizontal x-y average when the flow statistics are evolving in time, and with an additional time average over an inertial period when the flow is quasi-steady. It is important to note that u * and δ N are functions of the Reynolds number and are not known prior to simulation in Ekman layer studies. Denoted by superscript +, statistics normalized with the viscous scale (ν/u * ) are chosen to describe behaviour in the near-surface region. Normalization of statistics, denoted by superscript −, by the boundary-layer height (Ekman layer thickness δ N ) is also used.
2.1. Surface roughness The roughness takes the form of periodic two-dimensional periodic bumps whose non-dimensional amplitude is fixed and the steepness is changed among cases. The rough surface is described by a height function, η(x), which is generated through a harmonic function, f (x) = −h cos(2πx/λ), with wavelength λ and amplitude h, including only the positive portions of f (x), (2.5) A schematic of the roughness (figure 1) shows the surface has bumps but, unlike a surface wave, it has no troughs. This choice is motivated by the example of a two-dimensional hilly surface. Other types of surface roughness, e.g. a wavy surface that has been studied in the past, are worth future consideration. The chosen roughness has not only a small height (h + = 15) relative to the Ekman layer thickness but also has a gentle slope (h/λ 1), so that, in the unstratified cases, its effect is small. In the stratified cases, as will be shown and explained, the same roughness height has a strong effect that substantially changes the flow with respect to its smooth-bottom counterpart. The influence of roughness is found to extend up to the region (Ri g > 0.1) where buoyancy alters the mean flow and turbulent stresses.

Buoyancy
Buoyancy is imposed by a surface flux, whose strength is quantified by the Obukhov length, where κ is the von-Kármán constant and q 0 = −α∂ z θ | z=0 is the applied surface cooling flux. Here L provides an estimate of the height at which the buoyancy flux and turbulent energy production by mean shear are balanced. Using inner-layer scaling, the normalized Obukhov length becomes The gradient Richardson number, a function of z and t in this flow, is defined by where N 2 = gβ∂ θ /∂z is the squared buoyancy frequency and S 2 = (∂ u /∂z) 2 + (∂ v /∂z) 2 is the squared mean of the vertical shear. Large local values of Ri g imply suppression of shear production of turbulence, and Ri g = 0.25 is a stability boundary for the stratified shear layer. The non-dimensional inverse Obukhov length scale can also be interpreted in terms of Ri g at the surface, where Pr = 1 has been assumed. An alternative normalization of L is based on outer-layer coordinates, where δ N = u * /f is the boundary-layer scale. The neutral EBL has a boundary-layer height of approximately 0.5δ N . Since u * is a time-dependent function in the stratified DNS cases, so are L, L + and L − . The bulk Richardson number, an overall stability measure of the flow, is defined as where θ 0 = θ ∞ − θ 0 is the difference between the temperature above the boundary layer and the surface temperature. In the present study the surface temperature (θ 0 ) is a time-dependent variable during the time interval, T, over which finite-time constant-flux stability is imposed. The value of Ri b is also initially dependent on time but it becomes constant at the end of the time interval, T.

Numerical details
The governing equations (2.1) and (2.2) are numerically advanced in time using a combination of the low-storage third-order Runge-Kutta (RKW3) and mixed spectral-physical spatial discretization. The equations are written in generalized coordinates and the grid conforms to the bottom wall as described by Gayen & Sarkar (2011). Spatial discretization and derivative calculations in the spanwise direction are performed using Fourier transforms, and the derivatives in the streamwise and vertical directions are computed using a second-order central difference scheme. The nonlinear advection terms are dealiased with the 2/3 rule and a sharp-cutoff filter. The kinematic pressure, p, is computed by solving the Poisson equation that results from imposing zero velocity divergence at each time step. The boundary conditions for the velocity are no-slip and impermeability (u i = 0) at the surface, periodicity in the horizontal directions and stress free (∂u i /∂z = 0) at the upper boundary. The temperature gradient is fixed at the bottom surface to impose the desired cooling flux. A sponge region with Rayleigh damping is applied to u i and θ to minimize the spurious reflection of gravity waves in the upper boundary (z = L z ). The in-house solver, which has been developed for environmental flows, has been applied to the Ekman boundary layer (Gohari & Sarkar 2018) as well as complex geometries including stratified oscillating flow over a slope (Gayen & Sarkar 2011) and a triangular ridge (Rapaka, Gayen & Sarkar 2013;Jalali, Rapaka & Sarkar 2014).
A cooling flux, as determined by L + , is imposed in the neutral cases too in order to simulate the behaviour of a passive scalar. All of the rough and smooth-bottom stratified cases are initiated with a fully developed velocity field taken from the corresponding neutral case. The passive-scalar temperature field from the neutral case is reset to a uniform Case Here N x , N y and N z are the number of grid points in the streamwise, spanwise and vertical directions, respectively. The case label in column 1 has a subscript N (neutral case) or S (stable case), and starts with the number of bumps in the rough cases. In the stable cases, a surface buoyancy flux, chosen to obtain a target L + ≈ 700, is applied for a finite time of f T ≈ 6 and then the surface temperature is held constant.
background value so that each stratified case starts with zero temperature variation and stratification is allowed to build up in response to the applied surface buoyancy flux. Parameters used in this DNS study are summarized in table 1. A fixed value of buoyancy flux, corresponding to L + ≈ 700, is applied for a time interval of f T = 6 in the stable cases. Note that the computational domain is enlarged to twice the streamwise domain of Gohari & Sarkar (2018) in order to better accommodate long streamwise structures. The resolution is x + ≈ 8.5 in the streamwise direction for the flat-bottom case and y + ≈ 5.2 in the spanwise direction, similar to other DNS of wall-bounded flows. For the rough cases, the resolution in the streamwise direction is increased to x + ≈ 5.6. In the vertical direction, ten grid points span 0 < z + ≤ 10 with a non-dimensional grid spacing z + min ≤ 1. The roughness height is kept constant at a small value of h + = 15 which corresponds to the transitionally rough regime. The roughness takes the form of a periodic array of two-dimensional, spanwise-uniform elements whose surface elevation is given by (2.5). The roughness amplitude, h + = 15, is kept constant, and the wavelength of the roughness (λ) is changed. Note that λ = L x /N b , where N b is the number of bumps and L x is the streamwise domain size. In the simulations L x is kept fixed and λ is changed by changing N b . Doubling N b decreases the element half-length (l = λ/4) by a factor of 2 and doubles the slope. The coverage of the bottom by the roughness elements is 50 % of the wall area, independent of the number of bumps. The slope of each bump, given by h/l, is small and changes from 0.042 to 0.084 when N b is doubled from 2 to 4. Another measure is the maximum slope, hk where k = 2π/λ is the wavenumber, which changes from 0.06 to 0.12. We find that the flow does not separate at these values of slope. It is worth noting that the slope utilized here is below the critical slope for separation reported in the literature on a sinusoidal wavy surface where Gong et al. (1996), Sullivan et al. (2000) quote hk > 0.3 for flow separation, and Zilker et al. (1977) state that there is incipient separation at 2h/λ = 0.05 and there are large separated regions at 2h/λ = 0.125 and 0.20.

Flow evolution
The effect of the small-amplitude bumps on the velocity and temperature profiles in z is found to be insignificant under neutral conditions in contrast to the substantial effect under stratified conditions. This primary result is elaborated and explained in this section. Figure 2 shows the horizontal mean velocity ( u 2 + v 2 /U ∞ ) profile of the EBL in the top row and the normalized potential temperature (u * N (θ − θ ∞ )/q 0 ) in the bottom row for both neutral (left column) and stratified (right column) cases. The statistics are obtained by averaging over the horizontal x-y plane and an average over one inertial period ( ft ≈ 2π). The temperature is treated as a passive scalar in the neutral cases. The velocity in the neutral flat-bottom case compares well with previous results (Shingai & Kawamura 2004;Ansorge & Mellado 2014;Shah & Bou-Zeid 2014) at comparable Re as discussed by Gohari & Sarkar (2018). The horizontal wind speed exhibits little difference among the neutral (figure 2a) flat and rough cases. Since the slope of the bump is gentle (hk = 0.06 and 0.12), the flow does not separate and the change in wall drag (shear stress plus form drag) is small. It is worth noting that DNS of a turbulent flow over a sinusoidal wavy wall with hk = 0.1 in Couette flow (Sullivan et al. 2000) and channel flow (De Angelis et al. 1997) does not show flow separation. The stratified cases also do not exhibit flow separation; however, the mean velocity profiles (figure 2b) show a strong influence of roughness on the features of the LLJ that forms in the boundary layer. Each stratified case has a super-geostrophic velocity, commonly referred to as the LLJ. The formation of the LLJ is a distinctive feature of the stable ABL (Beare et al. 2006) associated with the reduction of turbulent momentum flux by buoyancy. In the roughness layer at the surface, the bumps counteract the buoyancy-induced reduction of the momentum fluxes. Thus, the peak of the LLJ moves upward and the LLJ profile broadens. In the stratified 4Bump case (filled diamonds in figure 2b), the wind-speed profile in the region between the surface and z − = 0.1 is close to the neutral case (open diamonds) while, in contrast, the 2Bump and flat cases exhibit a significant deviation from the neutral case. The present DNS results show a LLJ with a nose at z ≈ 0.1δ N , with a maximum super-geostrophic overshoot of u/u ∞ − 1 ≈ 10 %. Direct numerical simulation with stronger stability conducted by Gohari & Sarkar (2017) showed cases with LLJ at z ≈ 0.05δ N and u/u ∞ − 1 ≈ 50 %.  (2003) studied sound propagation in a stable nocturnal boundary layer which had a deep temperature inversion and LLJ at approximately 160 m, and was observed during CASES-99 (Poulos et al. 2002). It is worth noting that sloping terrain that leads to drainage flows also contributes to the LLJ structure (Mahrt 1999).

Overall velocity and temperature variation
In the neutral EBL (figure 2a) the effect of the surface roughness on the mean temperature is relatively small, similar to that on the velocity. However, in the stratified boundary layer (figure 2b) the bumps have a significant effect on the temperature distribution. The strong near-surface inversion of the flat case is substantially weakened in the 4Bump case; the near-surface temperature field is more mixed, and its profile moves towards that of the neutral case. Thus, in spite of employing the identical value of surface cooling flux (q 0 ) in the three stratified cases, the surface temperature in the 4Bump case does not decrease as much as in the other stratified cases.  The mean velocity is also examined using inner-layer coordinates. A least-squares fit to the vertical variation of G = ( u 2 + v 2 ) 1/2 is performed to obtain the following profile in semi-logarithmic coordinates: Here z 0 + = e −κB . The values of (κ, z 0 + ) are (0.44, 0.0759), (0.43, 0.0814) and (0.40, 0.1187) in the flat, 2Bump and 4Bump cases, respectively, and the profiles are as shown in figure 3. Sullivan et al. (2000) in their DNS of Couette flow found that, for a stationary bottom wall, z 0 + = 0.17 in the flat-bottom case and z 0 + = 0.27 for a wavy-bottom wall with ak = 0.1. Both cases exhibited κ = 0.41.
Monin-Obukhov (MO) similarity theory, often used to interpret ABL data, is utilized to assess the mean velocity and temperature profiles obtained here. Simulation data are used to compute stability functions (Φ m and Φ h ) associated with MO similarity theory: Here u * (z) and θ * (z) are the local scales for velocity and temperature fluctuations, as per local similarity theory (Nieuwstadt 1984): In the unstratified cases, Φ h is computed using passive-scalar statistics (buoyancy term in the momentum equation is set to zero) and the computed value of L L is to be understood as a notional value that allows comparison of profiles with the stratified cases. Figure 4 shows normalized mean gradients (the so-called stability functions) for unstratified (a,b) and stratified (c,d) Ekman layers. For the unstratified cases, Φ m (z) = 1 is expected in the log-law region and, correspondingly, Φ m (z) takes values near unity over an extended region (figure 4a). Very near the bottom and in the roughness sublayer, Φ m (z) increases with increasing z. According to MO theory, Φ m and Φ h are constant and close to unity when z/L L 1; here L L is the local Obukhov length. When z/L L ≥ O(1), the turbulent length scale becomes limited by the local Obukhov length and Φ m (z) increases with z. As shown in figure 4(c,d), there is a region, 0.02 < z/L L < 0.1, where Φ m is approximately constant but greater than unity, followed by an increase of Φ m as a function of increasing z/L L . The function Φ h (z) also increases with increasing z/L L and exhibits a slope that is larger than that of Φ m (z). Thus, the effect of buoyancy on heat transport is stronger than on momentum transport, i.e. the turbulent Prandtl number becomes larger than 1 in the stratified region of the boundary layer. The dependence of the stability functions on z/L L in the present work is similar to that reported in previous studies, e.g. the LES of Basu & Porté-Agel (2006). It is worth noting that, in the surface layer and among the stratified cases, the 4Bump case shows behaviour closest to the passive-scalar counterpart.   Figure 5 shows the overall influence of the bumps on the flow evolution. The integrated TKE, obtained by a horizontal x-y average to compute u i u i followed by integration in the vertical, is defined as In all cases, TKE initially decreases during a period of turbulence collapse when the flow transitions from neutral to stable. For the flat case, turbulence collapses with a time scale of L/u * in agreement with Flores & Riley (2011). The collapse is followed by a recovery of TKE in each case. The turbulence recovery is consistent with DNS results of the stable EBL (Ansorge & Mellado 2014;Shah & Bou-Zeid 2014;Gohari & Sarkar 2018). The value of u * decreases by ∼10-15 % during the initial transient before recovering to approximately its initial value. In an analysis of the CASES-99 data, Banta et al. (2007) reported a 6 hr collapse time period which corresponds to a non-dimensional time of ft ≈ 1.98, which is similar to the DNS collapse time scale.
The TKE evolution after collapse exhibits significant differences among cases. In figure 5(a,c) for the flat and 2Bump cases, TKE exhibits a large-amplitude inertial oscillation with period, ft ≈ 2π. For the 4Bump case (figure 5e), the periodic modulation of TKE is not observed, but TKE exhibits a gradual increase (on average) and seems to reach a plateau beyond ft ≈ 25. We will discuss reasons for the difference in TKE evolution later.
We apply the same buoyancy flux for all cases; however, the final Ri b (right column of figure 5) is different among the three cases. The final value of Ri b for the flat, 2Bump and 4Bump cases is 0.485, 0.379 and 0.312, respectively. Thus, the modification of the flow by the roughness elements is sufficiently strong in the 4Bumps case, despite the small bump height and the gentle slope of the bump, to significantly decrease Ri b . The decrease in Ri b suggests that the buoyancy effect in the 4Bump case on turbulence is weaker, as will be demonstrated by quantification of turbulent fluxes.
The contour of local gradient Richardson number, Ri g defined by (2.8), is a measure of the height-dependent strength of static stability relative to shear instability, and is depicted in figure 6. There is a region extending up from the wall which is subcritical (Ri g < 0.25). The height at which Ri g crosses the critical value of 0.25 increases with the number of bumps so that the subcritical region of the Ri g profile expands significantly. The subcritical region that starts at the bottom reaches up to z − ≈ 0.2 in the 4Bump case instead of z − ≈ 0.075 in the flat case. The implication is that roughness changes the stability of the near-bottom flow to make it more vulnerable to shear instability. Thus, the behaviour of both Ri b and Ri g (z) suggest that the roughness elements, albeit small, substantially mitigate the stabilizing effect of buoyancy. In figure 7 instantaneous vertical vorticity contours for the three stratified cases at different times ( ft ≈ π in the left column and ft ≈ 4π in the right column) are shown on a horizontal plane close to the wall (z + ≈ 16) and near the crest of the bumps. These times correspond to points (a-f) on the TKE profiles shown in the left column of figure 5 and also to panels (a-f ) in figure 7. Comparison of the points demarcated as b, d and f on the time histories in figure 5 show that, at ft ≈ 4π, the integrated TKE is the same for 2Bump and 4Bump cases and is slightly smaller in the flat case. However, on comparison of figures 7(b), 7(d) and 7( f ), we find that the near-wall structures are substantially different, reinforcing the fact that an overall statistical measure of turbulence does not necessarily reveal the full picture of the flow state. In particular, with an increasing number of bumps, near-wall turbulence is less patchy and more continuous at ft = 4π, corresponding to a state of continuous turbulence without global intermittency (Nieuwstadt 1984). This is true even at the earlier time of ft = π during the initial adjustment of the boundary layer to buoyancy when the TKE drops.

Boundary layer thickness
Previous studies have chosen different metrics to quantify the thickness of the stratified boundary layers, including but not limited to the height of the capping inversion layer (Melgarejo & Deardorff 1974;André & Mahrt 1982), the height at which the low-level jet velocity is maximum (Blackadar 1957;Shapiro & Fedorovich 2010;Van de Wiel et al. 2010), and the height at which turbulent stress reduces to some fraction of its surface value (Zilitinkevich 1972;Businger & Arya 1975;André & Mahrt 1982;Kosović & Curry 2000). Although each of these definitions has a suitable use, the one defined based on the location where turbulent stress vanishes is chosen here as an average measure of the interface between turbulent and non-turbulent layers. We define the height by locating the position (denoted by z p ) where the horizontal Reynolds shear stress is reduced to 5 % of u 2 * , and then linearly extrapolating to the location at which it would vanish if the stress profile was linear. Thus, a time-evolving value of the stratified boundary-layer thickness is defined as Subsequently, a modified bulk Richardson number, based on the local (in time) boundary-layer thickness, is defined as (3.8) , and 4Bump . The roughness element height is also shown in (a).
Figure 8(a) shows the time evolution of δ t for the stratified cases. In the flat case, the EBL thickness decreases sharply during the initial collapse of turbulence. Although small relative to the initial neutral boundary-layer height, the reach of the roughness bumps becomes comparable to the reduced value of δ t during turbulence collapse. Therefore, the roughness elements are able to sufficiently perturb the thin, collapsing boundary layer to partially arrest turbulence collapse. The modified bulk Richardson number (Ri b,t ) is shown in figure 8(b). The previously shown Ri b , based on the neutral boundary-layer scale (δ N ), is also shown for ease of comparison. For the 4Bump case, the non-steady values of Ri b,t are initially higher when the flow goes through turbulence collapse, however, the final values are similar. Thus, the overall strength of stratification as measured by Ri b,t does not change among cases. It is the wall-normal distribution of temperature and velocity which is affected by roughness. The invariance of Ri b,t among cases implies that δ t ∝ θ 0 −1 . Evidently, the introduction of surface bumps decreases the net amount of surface cooling (shown by the decrease of Ri b ) and, concurrently, increases the boundary-layer thickness to maintain Ri b,t . It is worth noting that, in previous DNS of the stable flat-bottom case with constant temperature boundary condition that imposes Ri b , the vertical extent of the Reynolds shear stress profiles also exhibits a similar trend of δ t increasing with decreasing Ri b . For example, it can be inferred from figure 13(c) of Shah & Bou-Zeid (2014), which shows the Reynolds shear stress for various stability levels, that δ t is approximately inversely proportional to Ri b .

Turbulent fluxes
Roughness enhances turbulent fluxes and, furthermore, the increase is substantially stronger in the stratified situation relative to its unstratified counterpart. Figure 9(a,b) shows profiles of the turbulent momentum flux ( u w 2 + v w 2 /u 2 * N ), which are obtained by horizontal x-y averaging and a time average over ft ≈ 2π. In the neutral EBL (figure 9a) the increase with respect to the flat case is negligible for the 2Bump case and moderate for the 4Bump case. The increase is substantially more in the stratified cases (figure 9b), e.g. the peak value in the 4Bump case is twice that in the flat case. Surface cooling suppresses the turbulent momentum flux as revealed by comparison of figures 9(a) with 9(b). However, surface roughness is able to counteract the suppression over a significant portion of the initial neutral boundary layer. The flux profiles in the stratified rough cases have a larger vertical extent than the stratified flat case, consistent with the increase of the boundary-layer thickness (δ t ) induced by the bumps. Buoyancy flux, w θ , in the stratified cases (figure 9d) is considerably suppressed relative to the corresponding unstratified cases (figure 9c), regardless of the number of bumps. After its peak, the magnitude of w θ drops sharply with increasing height relative to the unstratified cases. The drop is less sharp in the presence of roughness. Thus, between z − of 0.15 and 0.3 in figure 9(d), the value of w θ is substantially larger in the 4Bump case (filled red diamonds) relative to the flat case (filled black squares).
Roughness preferentially enhances vertical fluctuations in the stratified EBL. Figure 10 depicts the effect of roughness on the amplitude of horizontal (u h,rms = u 2 rms + v 2 rms ) and vertical (w rms = √ w w ) fluctuations. The presence of surface roughness leads to a mild increase of both u h,rms and w rms above the bumps in the neutral cases (left column). Under stratification, the surface roughness effect on w rms is dramatic. In the near-surface region (z − < 0.15) the 4Bump case has substantially larger w rms relative to the flat case, as seen in figure 10(d). This increase of vertical transport is key to the roughness-induced increase of TKE, as will be evident later in § 3.4 where the TKE budget is discussed. It is worth noting that the increase of near surface w rms in the 4Bump stratified case is also manifested in the finding, illustrated by figure 7( f ), that the 4Bump case has continuous near-surface turbulence in contrast to the local intermittency (localized turbulence patches distributed in an almost-quiescent background) of the flat case in figure 7(b). The location of the EBL boundary (z = δ t ), based on the turbulent momentum flux, is shown on each profile in figure 10. The r.m.s. fluctuation profiles have significantly larger vertical spread than δ t which is based on the Reynolds shear stress. Mahrt (1998), based on observational data, presents idealizations of the stratified boundary layer that we assess in figure 11 using the present simulation data. The very stable boundary layer (the right subfigure of figure 1 in Mahrt 1998) is conceptualized by the author as follows: a w rms profile that has weak near-surface values with the peak occurring at an elevated location, and a θ profile that exhibits a strong surface inversion layer. The flat case (figure 11b) shows a strong near-surface inversion and an elevated peak of w rms (associated with the shear of the LLJ), in good agreement with the idealization of a very stable boundary layer. On the other hand, the 4Bump stratified case (figure 11a) has a presentation that is more akin to the so-called weakly stable boundary layer (the left

Turbulent kinetic energy budget
The TKE budget is analysed to better understand the mechanisms underlying turbulence collapse and rebirth. At statistical steady state, the Reynolds-averaged turbulent kinetic energy budget can be written as where turbulent kinetic energy e = 1 2 u i u i , pressure transport rate P T = −∂ u i p /∂ x i , turbulent transport T T = − 1 2 ∂ u i u i u j /∂ x j , shear production rate P = − u i u j ∂ u i /∂ x j , buoyancy flux B = δ i3 βg u i θ , viscous diffusion rate ν D = ν∂ 2 e/∂ x 2 j and viscous dissipation rate = ν ∂u i /∂ x j ∂u i /∂ x j . It is worth noting that the buoyancy flux (B) is negligible in comparison to the other terms in the TKE balance. The smallness of B in the stratified boundary-layer TKE budget is consistent with other studies of the EBL, e.g. Shah & Bou-Zeid (2014) and Gohari & Sarkar (2018) who found that B is negligible in the TKE balance and that the direct impact of stability on the flow is not through TKE destruction by buoyancy, but rather through the inhibition of shear production. The balance, calculated as the sum of all the terms on the right-hand side of (3.9), has been obtained as a function of time. In the 4Bump case the sum is less than 1 % of the dominant term. In the flat and, to a lesser extent in the 2Bump case, the instantaneous sum is oscillatory owing to the inertial oscillations in TKE. The amplitude of these oscillations can be as large as 15 % but the time average over an inertial period is again less than 1 % of the dominant term. Turbulent kinetic energy budget terms are shown in figure 12 for flat (a-c), 2Bump (d-f ) and 4Bump (g-i) cases at ft = 0 (a,d,g), ft ≈ π (b,e,h) during the initial turbulence collapse, and ft ≈ 3π (c, f,i) after turbulence recovery. At ft = 0, which corresponds to the neutral EBL, the TKE budget terms are similar among the three cases. The peak of TKE production (P) is in the buffer layer. After the imposition of constant-flux stability, there is a drop of P. For the flat case (figure 12b), P becomes negligible at ft ≈ π. The restriction of TKE production in the surface layer by the imposed stable surface cooling flux leads to turbulence collapse (Ansorge & Mellado 2014;Shah & Bou-Zeid 2014). In the flat case, Gohari & Sarkar (2017) showed that turbulence recovery is promoted by pressure transport (P T ) that carries outer-layer fluctuation energy into the lower flank (with subcritical Ri g ) of the near-surface LLJ. The pressure transport is important during the recovery in the flat case (green line in figure 12c), but not so in the cases with roughness ( figure 12f,i). For the 2Bump case in figure 12(e), the magnitude of the shear production is reduced, but nevertheless, roughness ensures that |P| > 0. The production of mechanical turbulence enhances the subsequent recovery of the EBL to a turbulent state. For the 4Bump case in figure 12(h), the reduction of shear production by buoyancy is even weaker than the 2Bump case. It is evident from figure 12 that the surface roughness aids TKE generation to keep the flow from turbulence collapse. For the same reason, near-surface turbulence is continuous in the 4Bump case after recovery rather than being locally intermittent as in the flat case. The reduction of P by buoyancy is related to the significant damping of vertical turbulent motions which in turn reduces the momentum fluxes ( u w and v w ), especially in the flat case, as seen by comparing figure 9(b) with figure 9(a). This reduction is mitigated in the rough cases because the surface bumps promote vertical transport.

Coherent structures
The roughness elements introduce coherence into the flow vorticity and velocity. Figure 13 shows a λ 2 isosurface superposed on the streamwise vorticity (ω x ) contour on a near-surface horizontal plane. The quantity, λ 2 , is the intermediate eigenvalue of the symmetric tensor, S ik S kj + Ω ik Ω kj , where S ij and Ω ij are the symmetric and antisymmetric components of the velocity gradient tensor, ∂u i /∂ x j . The unstratified flow (left column) exhibits coherent packets of hairpin vortices. The near-wall structures are inclined with a veering angle with respect to the outer geostrophic velocity which points in the x-direction. The signature of the bumps is seen in the coherent strips of ω x on the displayed horizontal plane. The snapshots are shown at ft ≈ 6 when the TKE in the stratified cases is approximately at its minimum. At this time, the coherent structures in the stratified EBL (right column of figure 13) are suppressed relative to neutral conditions. Nevertheless, the roughness elements are able to sustain some of the unsteady flow structures as can be seen by comparing the 4Bump case (figure 13f ) with the flat case (figure 13b).

Dispersive effects of roughness
The roughness elements introduces a variability in the flow relative to the horizontal average. This variability leads to a so-called dispersive component of velocity which brings about fluxes of momentum and heat, and is extracted as follows. First, the field variables are decomposed into a time-averaged mean and a fluctuating component. Thus, the velocity on a given horizontal plane (fixed z) is split into whereū i is the time-averaged mean and u i is the fluctuation about the time-averaged mean. The two-dimensional (2-D) periodic roughness in the present problem imposes a spatial organization on the time-averaged field that can be extracted by applying a further streamwise spatial average, denoted by x , to the time average (ū i ). Thus, (4.1) becomes This decomposition identifiesũ i (x, y) as the spatially coherent part of the time-averaged flow, e.g. Finnigan (2000), Sullivan et al. (2000), Poggi et al. (2004), Li & Bou-Zeid (2019). Physically,ũ i (x, y) in the present problem is the time-mean velocity associated with the roughness bumps. Because of the streamwise periodicity and statistical steadiness, the composite (t and x) average, ū i x ( y), is taken to be the Reynolds average Isosurface of λ 2 = −3.125(u * N /z) 2 is shown. and the fluctuation, u i , is the Reynolds fluctuation. As elaborated in this section, the u i field contains a time-independent part (ũ i (x, y)) associated with the 2-D roughness elements which leads to a dispersive component of the Reynolds stress and an 'incoherent' part (u i (x, y, t)). The triple decomposition (Reynolds & Hussain 1972) splits the velocity into a time-average, a time-coherent (often a wave-like field at a specific temporal frequency) component, and an incoherent turbulent component. The first line of (4.2) is analogous to their triple decomposition as it extracts the space-coherent component introduced by the periodic roughness.
For completeness, the mathematical definition of the averages and the fluctuations are given below for the velocity on a horizontal plane at constant z: (4.7) Figure 14 shows the consequences of the decomposition of streamwise velocity (u) in the neutral EBL. The neutral value (u * N ) of friction velocity is used to normalize velocity in figures 14-18. Upon comparison of figure 14(a,c), it is evident that the Reynolds fluctuations (u ) and the incoherent velocity fluctuations (u ) are similar in the unstratified flat case andũ = 0. However, in the 4Bump case the u field (figure 14b) is different from the incoherent u field (figure 14d). The u field contains a spatially organized component, which is introduced by the bumps. This coherent component (ũ i ), isolated by applying the triple decomposition, appears as four distinctive strips associated with the four surface bumps in figure 14( f ). The region of positive u which is forward of a surface bump combines with the rearward negative u to form a strip associated with the roughness. The imprint of the 2-D bumps, extracted throughũ, is unequivocal.
The influence of the bumps extends upward from the roughness elements into the flow as illustrated by figure 15. The vertical coherent component (w in figure 15c) takes values of O(u * ) up to z − = 0.1, which is ∼3 times the roughness height. When the boundary-layer thickness (δ t ) decreases during the initial transient as the EBL adjusts to the imposed cooling flux, the vertical reach of the bump-induced flow is no longer small compared to δ t . For example, during the initial turbulence collapse in the flat case, the boundary layer thins to δ − t ≈ 0.1 (figure 8a). The vertical extent of the roughness-induced perturbation to the flow is sufficient to mitigate the turbulence collapse and δ − t does not become smaller than 0.2 in the 4Bump case.
The dispersive effect of the roughness is an important contributor to the Reynolds shear stress as is demonstrated by figure 16. Let M = u w 2 + v w 2 denote the turbulent momentum flux obtained after Reynolds averaging. Since the Reynolds fluctuation can be decomposed as u =ũ + u , it follows that (4.8) where the coherent (M coh ) and incoherent (M inc ) components are given by (4.9a,b) and the cross-term is (4.10) In figure 16(a), 4Bump has a higher peak value of M relative to the flat case. The incoherent part (figure 16b) has similar values for the peak, independent of surface roughness. The difference in M among cases arises from the coherent part (figure 16c) and also the cross-term (figure 16d). Both, coherent and cross-terms, are negligible in the flat case, and their values in the surface layer increase with an increasing number of bumps. The preceding discussion shows that the presence of surface bumps substantially modifies the instantaneous velocity and the turbulent momentum fluxes by introducing a dispersive component. In the stratified cases the dispersive component remains substantial and, furthermore, its importance to the turbulence energetics is increased since the incoherent part is suppressed by buoyancy. The instantaneous Reynolds shear stress (u w ) carries contributions from (ũ,w) alone, (u , w ) alone, and their cross-correlations, as illustrated for the 4Bump stratified case. The initial turbulence collapse ends when ft ≈ 3 has elapsed after imposing the surface cooling flux. The instantaneous Reynolds shear stress (u w in figure 17a) at ft ≈ 3 shows the presence of the 2-D coherent bumps which is comingled with the incoherent fluctuation (u w in figure 17b). The dispersive component (u w − u w in figure 17d) is substantial.
In the stratified cases roughness enhances the turbulent fluxes relative to the flat case by not only introducing coherent (ũw) and dispersive (u w − u w ) components into the shear stress but also by significantly changing the incoherent component (u w ) relative to the flat case. Figure 18 shows the Reynolds shear stress (same as the incoherent component in the flat case) at early and late times in the flat, stratified case. Direct comparison of figure 17(b) with figure 18(a) shows that, at the early time of ft ≈ 3, the 4Bump case has significant near-wall Reynolds shear stress in contrast to the negligible level in the flat case. Later in time, u w recovers somewhat in the flat case (figure 18b) although it is spatially sparse relative to the rough case.

Summary and conclusions
The stratified Ekman boundary layer has served as a canonical problem that is amenable to high-resolution simulation and is relevant to the stratified ABL. However, roughness effects in the EBL have not been studied previously using DNS. This motivates the present investigation of an EBL on a surface with two-dimensional bumps. A DNS study is conducted where Re * is fixed at approximately 700, the cooling flux is fixed at a non-dimensional value of L + ≈ 700 which is sufficient to induce turbulence collapse during the initial transient (Gohari & Sarkar 2018), the roughness amplitude is fixed at a small value in the transitionally rough regime, and the slope of the elements (equivalently, number of bumps per unit length) is varied. The focus is on the competition between flow stabilization by buoyancy and possible destabilization by the roughness.
We find that the flow evolution is substantially affected by roughness in the stratified EBL, especially so in the case with 4Bumps, the case with the highest number of bumps per unit length and also the highest geometrical slope. In the 4Bump case the minimum TKE reached during the initial turbulence collapse is significantly larger than in the flat case. Furthermore, later in time after turbulence recovery, the near-surface flow state is continuously turbulent in the 4Bump case in contrast to the local intermittency (turbulent patches interspersed within quiescent regions of near-laminar flow) in the flat and 2Bump cases. The MO stability functions show that, at locations where z is not small compared to the local Obukhov length, buoyancy affects the mean momentum and temperature profiles. Furthermore, this buoyancy effect is substantially weaker in the 4Bump case although the applied surface cooling flux is identical among cases. By increasing the number of bumps per unit length keeping bump height constant, the slope of the roughness element is systematically increased in the present DNS. It is the increase in slope that enhances the magnitude of the roughness-associated effect on turbulence from 2Bump to 4Bump cases. Mahrt (1998) classifies the stable ABL into two regimes: (i) a very stable ABL that has a strong thermal inversion at the surface and vertical fluctuations which, although very small near the surface, have a peak at an elevated location; and (ii) a weakly stable ABL that has a weak thermal inversion and vertical fluctuations that, although somewhat reduced near the surface with respect to the neutral state, are substantial and do not display a peak at an elevated location. Examination of mean and turbulence profiles in the DNS results reveal that the flat case belongs to the very stable regime while the 4Bump case belongs to the weakly stable regime. Mahrt (1998) distinguishes between these two types of stable boundary layer based on environmental forcing: the weakly stable case is associated with windy night-time conditions, and the very stable case is associated with calm conditions and clear night-time sky. Here, we find that, for the same environmental forcing, roughness elements induce a qualitative change in the boundary layer, namely, from the strongly stable regime to the weakly stable regime. As the neutral BL thins and near-surface turbulence weakens during the initial response to the applied cooling flux, the reach of the roughness elements is sufficient to keep the local gradient Richardson number from becoming supercritical (exceeding 1/4) and concurrently maintain turbulence structures. The present result is obtained for a relatively low-Re flow amenable to DNS, where the roughness elements correspond to a transitionally rough regime. It is possible that, at higher Re and in a fully rough ABL, roughness elements can also lead to a qualitatively similar destabilization if their reach becomes sufficient compared to the buoyancy-induced thinning of the ABL.
The mechanism of turbulence collapse and rebirth is investigated in detail with the analysis of turbulent fluxes. The inhibition of surface-layer shear production, associated with the buoyancy-induced reduction of the turbulent momentum flux, is the reason for turbulence collapse as was found by Ansorge & Mellado (2014) ;Shah & Bou-Zeid (2014). The positive buoyancy flux in the TKE budget is not the reason. Collapse is followed by turbulence recovery in both flat and rough cases. Notably, the mechanism of turbulence recovery in the rough cases is different from that in the flat case. In the flat case a strong LLJ is formed at the end of turbulence collapse and, as found by Gohari & Sarkar (2017), pressure transport of fluctuations from the outer layer into the sheared lower flank of the LLJ triggers locally intermittent turbulence as the fluctuations interact with the enhanced shear of the LLJ. In contrast, we find in the 4Bump stratified case that the enhanced vertical (w) velocity in the near-surface roughness layer counteracts the buoyancy-induced reduction of turbulent momentum flux. Unlike the flat case, shear production of TKE does not decay to a negligible value in the rough cases during the initial collapse of turbulence.
The roughness elements introduce a spatial organization into the flow. A triple decomposition is employed to isolate the coherent component. The roughness-associated dispersive component (total minus the incoherent part) of the Reynolds stress is found to be substantial in both unstratified and stratified cases. Since the incoherent component is strongly suppressed by buoyancy, the dispersive component becomes more important to the turbulence dynamics in the stratified cases. Additionally, in the stratified cases, the incoherent component is enhanced with respect to the flat case.
The present small-amplitude bumps of h + = 15 correspond to a transitionally rough regime, have a gentle slope which does not lead to flow separation, and have a small effect on the flow in the neutral EBL. However, since the layer of boundary-associated turbulence thins during the initial collapse, the bump height becomes sufficient for the influence of roughness to reach into an appreciable portion of the boundary layer so as to modify the shear and stratification in the boundary layer. In particular, the near-bottom region of subcritical Ri g becomes substantially thicker in the rough cases and, correspondingly, the buoyancy-induced suppression of turbulent fluxes in the surface layer is mitigated. Thus, roughness modifies the boundary layer from a very stable regime in the flat case to a stable regime. For a sufficiently large cooling flux, larger than the value considered here, it is possible that the rough boundary layer reverts back to a very stable regime. In future work we will systematically vary the cooling heat flux and roughness element height to further understand the role of roughness in counteracting the effect of stable stratification. It will also be desirable to extend simulations in future work to higher Re and a fully rough regime. Another future direction is the consideration of more complex geometry, e.g. three-dimensional obstacles and multiscale roughness.

Declaration of interests
The authors report no conflict of interest.