Direct numerical simulation of turbulence over two-dimensional waves

Understanding turbulent boundary layers (TBLs) over complex surface undulations has immense practical value in both environmental and engineering applications. Common examples in engineering include internal flows in pipes and turbomachinery, external flows over fouled ship hulls,¹ and wind turbine blades. In the environment, the Earth’s surface offers a wide variety of surface heterogeneities including water waves, sand dunes, and shrubs (small); tree canopies and buildings (medium); and mountains and hills (larger scale). Understanding this ubiquitous flow dynamics over such a wide range of heterogeneities becomes a necessity. Given the interest in accurate drag prediction, a significant amount of research has been devoted to understanding turbulent flows over pipe roughness, for example, the work of Darcy² nearly two hundred years ago; in the early half of the 20th century by Nikuradse,³ Colebrook,⁴ and Moody;⁵ and more recently by various research groups.^6–8 By roughness, one typically refers to undulations (of size k) much smaller [say, k/δ ≲ O(0.01)] than the outer layer scale (δ, boundary layer thickness) and much larger [say, k/$Lv$ = k⁺ ≳ O(100)] than the inner layer scale (viscous length, $Lv$⁠). Turbulent flows over uniform roughness embedded in a flat surface have been studied extensively using laboratory experiments^9–16 and reviewed in Refs. 14 and 17. In addition, there exist studies of flow over systematically designed roughness using experiments⁹ and more recently using direct numerical simulation (DNS)^8,18–28 and large eddy simulation (LES).²⁹ These studies invariably focus on a relatively large blockage [k/δ ≲ O(0.01–0.1)] at transitional Reynolds numbers [k/$Lv$ = k⁺ ≳ O(10)] with a sufficiently large surface layer (low to moderate Reynolds numbers, Re_τ ∼ 100–700) to keep the computational cost manageable while sufficiently resolving the relevant flow features. Advances in mainstream computational algorithms have enabled interesting recent efforts to mimic Nikuradse-type sand grain roughness using DNS at moderately high Reynolds numbers^30,31 (Re_τ ∼ 400–700), making it a step closer to simulating realistic surfaces. Such capabilities make it possible to evaluate different surface textures using super computers to potentially identify optimal designs and also generate the next generation of industrial look-up tables for drag estimation.³² Complementing such efforts are simulation studies of turbulent flows over structured undulations that offer a more easily interpretable framework to relate statistical trends with physical mechanisms. Such studies help us to (i) characterize the roughness sublayer scale and its relationship to other flow scales, (ii) characterize the roughness layer dynamics for a given parameterized shape, (iii) characterize the interaction of sublayer dynamics with inertial layer turbulence, and (iv) develop strategies for physics-aware predictive models. Knowledge of the surface layer dynamics from such fundamental studies also serves the additional purpose of understanding large-scale surface heterogeneity-driven lower atmosphere dynamics such as those pursued by our group^33–35 and elsewhere.³⁶ It is in this context that the current work focusing on in-depth investigation of turbulence generation over parameterized two-dimensional wavy surfaces is undertaken. Therefore, insight from this work has implications both to roughness characterization over complex surfaces and to understanding surface layer dynamics in geophysical settings over terrain.

Early efforts^3–5 related drag with roughness that is broadly classified as hydraulically smooth, transitional, or fully rough (depending on the roughness scale) and corresponds to flows with dominant viscous, combined viscous and form drag, and dominant form drag, respectively. Therefore, in the fully rough limit with minimal viscous effects, the drag depends only on the roughness scale, k, and not on the Reynolds number. The classical view is that roughness impacts only the near surface turbulence structure (up to a few roughness lengths), resulting in the outer layer flow similarity or, more commonly, the Townsend wall similarity.^37,38 The roughness function,³⁹ Δ⟨u⟩⁺, quantifies the roughness induced downward shift in the outer layer of the mean velocity profile (when plotted in a semi-log scale) due to increased drag from surface heterogeneity. Topology impacts drag (and Δ⟨u⟩⁺) and also the surface layer dynamics. For example, two-dimensional wavy undulations generate stronger vertical disturbances (and Δ⟨u⟩⁺) due to the absence of significant spanwise motions in the roughness sublayer as compared to three-dimensional wavy surfaces.^20,22,40,41 Therefore, understanding how surface undulations impact the mean flow and consequently parameters such as Δ⟨u⟩⁺ is of great importance.

The current study is one step along this direction where we try to decipher the turbulence structure and production mechanisms over two-dimensional wavy surfaces. Turbulent flows over smooth sinusoidal surface undulations have been explored experimentally^42–47 since the 1970s with applications ranging from flows over water waves and dunes to roughness characterization. Zilker and co-authors^43,44 studied small amplitude effects of wavy surface abutting a turbulent flow, and note that the dynamics is well approximated by linear theory for cases with very little to no separation. Larger amplitude wavy surface turbulence⁴⁸ shows significant flow separation with strong nonlinear dynamics, especially for wave slopes, ζ = 2a/λ ≳ 0.05, where a is the wave amplitude and λ is the wavelength. This slope factor, ζ, can be generalized into solidity⁴⁹ for complex surface shapes. For flows with incipient separation or when fully attached (ζ = 2a/λ ≲ 0.03), the pressure variations show a linear dependence with the wave height. With advances in computing, DNS of turbulent flow over two-dimensional sinusoidal waves^18,20 has been increasingly leveraged to analyze the detailed turbulence structure at moderate Reynolds numbers. De Angelis, Lombardi, and Banerjee¹⁸ explored the asymmetric pressure field along the wave and analyzed the turbulence energy budgets. The near-surface dynamics is summarized as follows: The wavy surface undulations generate alternating and asymmetric bands of favorable (upslope) and adverse (downslope) pressure gradients, which, in turn, cause regions of alternating shear (impacting dissipation) and Reynolds stresses (impacting production) that complement each other. In the downslope, the higher momentum fluid goes away from the surface (ejection-like events), while in the upslope, it moves toward the wall (sweep-like events). This results in shape-induced turbulent mixing, increased dispersive (and Reynolds) stresses, and pressure–strain generation of spanwise turbulent motions (through splat events). In studying TBL over three-dimensional (3D) egg-carton shaped sinusoidal undulations, Bhaganagar, Kim, and Coleman²⁰ showed that surface undulations can lead to turbulence generation even for linearized flow dynamics. Wavy surface undulations impact the near-wall coherent structures in a manner consistent with the horizontal scale of the undulations as estimated using two-point correlation measures.^18,20 To our knowledge, most literature on turbulence analysis over wavy walls (or other shapes) focus on the double-averaged turbulence structure. A subset of these characterizes the streamwise varying (2D) statistical structure. An even smaller fraction explores the 2D structure of production, dissipation, and Reynolds stress transport. Consequently, there is little knowledge of the component level production mechanisms and their spatial structure that could enhance our understanding of drag generation mechanisms in such flows.

In this study, we explore the near-wall turbulence structure over two-dimensional wavy surfaces using direct numerical simulation of wavy channel flow at a friction Reynolds number, Re_τ = δ⁺ ≈ 180, with higher-order spectral-like accuracy⁵⁰ and surface representation using an immersed boundary method (IBM).^51,52 The focus of our current analysis is to better understand component level turbulence generation at low enough slopes where viscous and form drag are both important. For the sinusoidal two-dimensional surfaces considered in this study, the wave steepness or solidity, ζ, is varied over 0–0.044 while keeping the inner-scaled wave height, a⁺, constant. This ζ range represents the transition from attached flow to incipient to weakly separated flow. The rest of this article is organized as follows: In Sec. II, we describe the numerical methods, simulation design, quantification of statistical convergence, and validation efforts. In Sec. III, we present the streamwise averaged first order turbulence structure. In Sec. IV, we characterize the second-order turbulence structure, namely, the components of the Reynolds stress tensor and their generation mechanisms as modulated by the wavy surface. In Sec. V, we summarize the major findings from this study.

In this study, we adopt a customized version of the Incompact3D⁵⁰ code framework to solve the incompressible Navier–Stokes equation for Newtonian flow described in a Cartesian coordinate system with x, y, and z corresponding to streamwise, vertical, and spanwise directions, respectively. The skew-symmetric vector form of the equations are given by

and

where p represents the pressure field responding to the flow and Π the driving pressure gradient. For this incompressible fluid, the fluid density (ρ) is taken as one. The above system of equations are then rewritten to generate a separate Poisson equation for pressure.

The system of equations are advanced in time using a 3rd order Adams–Bashforth (AB3) scheme with pressure–velocity coupling using a fractional step method.⁵³ The velocity is staggered by half a cell to the pressure variable for exact mass conservation. A 6th Order Central Compact Scheme (6OCCS) with quasi-spectral accuracy is used to calculate the first and second derivative terms in the transport equation. The pressure Poisson equation (PPE) is efficiently solved using a spectral approach. The right-hand side of the PPE is computed using a quasi-spectral accuracy using 6OCCS and then transformed to Fourier space. To account for the discrepancy between the spectrally accurate derivative for the pressure gradient and a quasi-spectral accuracy for the divergence term, the algorithm uses a modified wavenumber.

A major downside of higher order schemes is the representation of complex geometry. In particular, the boundary conditions for higher order methods are hard to implement without loss of accuracy near the surface. In this work, we adopt an immersed boundary method (IBM) framework where the solid object is represented through a force field in the governing equations to be solved on a Cartesian grid. In this study, we leverage a higher order, direct forcing approach requiring reconstruction of the velocity field inside the solid object so as to enforce zero velocity at the interface. Therefore, this velocity reconstruction step is the key to success of this approach. The numerous different IBM implementations⁵² differ in the details of this reconstruction. In this study, we adopt the one-dimensional higher order polynomial reconstruction of Gautier, Laizet, and Lamballais⁵⁴ and refer to the work of Khan and Jayaraman³⁵ for a more detailed presentation of the method. The reconstructed velocity field is directly used to estimate the derivatives in the advection and diffusion terms of the transport equation, which is shown to be reasonably accurate as per Sec. II C. Although the IBM is the popular option to resolve complex surface undulations, one may alternatively use curvilinear grids²⁶ that offer better resolution of the surface. Our choice in this work is motivated by a couple of reasons. The first is numerical, i.e., the need for higher order advection and diffusion schemes that are harder to implement in curvilinear coordinates. The second is practical, namely, the need to extend these research studies to a more complex topology for a detailed comparison for which it is beneficial to keep the numerical algorithm invariant.

We carry out six different simulations of turbulent channel flows with flat and wavy surfaces of different steepness levels (ζ), but constant wave amplitude (a), as shown in Fig. 1. Figure 1(a) shows the entire channel with representative wavy surfaces that mirror each other. In Fig. 1(b), we show the close-up of the near-surface grid with the 2D immersed surface resolved by a grid of dimension 216 × 24. Here, ζ is the average steepness for this sinusoidal shape given by ζ = 2a/λ, where λ is the wavelength. In this study, ζ is varied from 0 to 0.044 with 0 corresponding to a flat surface and ζ = 0.044 corresponding to four sinusoidal waves over the streamwise length of the simulation domain. In all these cases, care was taken to ensure that the realized friction Reynolds number is sufficiently close to the targeted value of ∼180 by modifying the corresponding bulk Reynolds number. The bulk Reynolds number, $Reb=ubδν$⁠, for the flat channel is chosen as ∼2800. For the different wavy channel flows with the same effective flow volume and mean channel heights, using the same flow rate (or bulk Reynolds number) increases the friction Reynolds number, $Reτ=uτδν$⁠, due to the increase in u_τ with wave steepness, ζ. Therefore, to achieve a constant value of the friction Reynolds number, the bulk Reynolds number was appropriately reduced through an iterative process so that the increment in u_τ does not significantly affect the friction Reynolds number, $Reτ=uτδν$⁠. The simulation parameters for the different cases are summarized in Table I. Although one could perform these studies at much higher Reynolds numbers, the choice of Re_τ = 180 was chosen to balance computational cost and storage requirements, yet maximize resolution in the roughness sublayer. We anticipate the Reynolds number to primarily affect the extent of separation dynamics, namely, thinning of the shear layer at higher Re_τ and implications thereof. Although we observe this in our work, a detailed analysis is saved for future.

The simulation domain is chosen as 4πδ × 2.2δ × 4πδ/3 (including the buffer zone for the IBM), where δ is the boundary layer height set to unity for these runs. This domain size is much larger than that used by Jelly and Busse⁵⁵ and comparable to the full span design of MacDonald et al.²⁵ This volume is discretized using a resolution of 256 × 257 × 168 grid points, which is more than adequate for the purposes of this study. In the streamwise and spanwise directions, periodic boundary conditions are enforced while a uniform grid distribution is adopted. In the wall normal direction, no slip condition representing the presence of the solid wall causes inhomogeneity. To capture the viscous layers accurately, a stretched grid is used. The grid stretching in the inhomogeneous direction is carefully chosen using a mapping function that operates well with the spectral solver for the pressure Poisson equation. The different inner scaled grid spacings are also included in Table I. These data show that our near wall grid spacing may be slightly coarser in the vertical direction while using a smaller growth rate away from the wall as compared to the DNS of MacDonald et al.²⁶ We quantify the normalized vertical grid spacing with respect to the inner-layer length scale in Fig. 2(a) and the Kolmogorov scale, η, in Fig. 2(b) as a function of y⁺. We see that the near wall vertical grid spacing resolves the Kolmogorov scale, η, through the entire boundary layer. This grid design limits stretching near the surface and packs more points in the buffer layer without increasing grid cost. As we will note in Sec. III C, this choice is also justified by thickening of the buffer layer, especially over two-dimensional undulations. Given the higher numerical accuracy of our method compared to the DNS algorithms of MacDonald et al.²⁶ and Jelly and Busse⁵⁵ who used 2nd order schemes as compared to spectral accurate 6th order schemes in our work, the slightly coarse near wall grid distribution is not expected to impact the overall accuracy, especially when viewed in the context of the effective grid that accounts for numerical diffusion errors.

To quantify the statistical stationarity of the turbulence simulation data, we consider the streamwise component of the inner scaled mean (spatial and temporally averaged) horizontal stress (including the viscous and Reynolds stress), $τH,x+=⟨∂u∂y⟩x,z,t+−⟨u′v′⟩x,z,t+$⁠. Here, ⟨⟩_x,z,t represents the averaging operation with subscripts denoting averaging directions. In the limit of statistically stationary and horizontally homogeneous turbulence, $τH,x+(y)$ converges to a linear profile, $1−yδ$⁠, as derived from mean momentum conservation. Using this, we estimate a residual convergence error, $ϵRes=⟨∂u∂y⟩x,z,t+−⟨u′v′⟩x,z,t+−(1−yδ)$⁠, whose variation with y/δ is shown in Fig. 3. As expected, this error is sufficiently small for the flat channel (ζ = 0) with magnitudes approaching O(0.001–0.01) through the turbulent boundary layer (TBL). The plot also shows a similar quantification for wavy channel turbulence data with large residual errors near the surface indicative of deviations from equilibrium due to local homogeneity. Such quantifications also provide insight into the height of the roughness sublayer beyond which the mean horizontal stress approaches equilibrium values.

We perform a baseline assessment of computational accuracy for the turbulent channel flow using an immersed flat channel surface before adopting it for more complex surface shapes. We compare statistics from the current DNS with the well-known work of Kim, Moin, and Moser⁵⁶ regenerated by Moser, Kim, and Mansour⁵⁷ (MKM99 here onwards), which corresponds to a bulk Reynolds number, Re_b ≈ 2800, a mean centerline velocity Reynolds number, Re_cl ≈ 3300, and a friction Reynolds number, Re_τ ≈ 180. MKM99 used nearly 4 × 10⁶ (128 × 129 × 128) grid points and solved the flow equations by advancing modified variables, namely, wall-normal vorticity and Laplacian of the wall-normal velocity without explicitly considering pressure. They adopt a Chebyshev-tau scheme in the wall-normal direction, Fourier representation in the horizontal direction, and the Crank–Nicholson scheme for time integration. In our work, we adopt spectrally accurate 6th order compact scheme in space and a third order Adams–Bashforth time integration.^35,50 Figure 4 shows that the inner-scaled mean [Fig. 4(a)] and root mean square of the fluctuations [Fig. 4(b)] from the current DNS match those of MKM99.

The primary goal of this study is to explore the non-equilibrium, near-surface turbulence structure over systematically varied sinusoidal undulations with particular emphasis on delineating the shape dependent turbulence generation mechanisms. Naturally, this involves characterization of deviations in the (streamwise-averaged) first order turbulence structure from equilibrium phenomenology as evidenced in flat channel turbulence, assesses the extent of outer layer similarity, and tries to relate these observations to roughness induced turbulence effects as appropriate. For ζ ≪ 1 considered here, separation is minimal as shown in Fig. 5 using isosurfaces of instantaneous negative velocity. We note that separation is inconsistent but becomes prominent at higher ζ. The streamwise-averaged or “double-averaged” turbulence structure denoted by ⟨⟩_x,z,t is a function of solely the wall normal distance and refers to averaging along both homogeneous (z, t) and inhomogeneous (x) directions. Temporal (t) averaging is performed using 2500 three-dimensional snapshots over 20 flow through times for the chosen friction Reynolds number. Data are output every 200 time steps, which corresponds to nearly a snapshot of data every viscous time scale $ν/uτ2$ and two snapshots for every near-wall Kolmogorov time scale, $τη=ν/ε$ [see Fig. 2(c) for plots of τ$v$_,η/Δt ratios with y⁺], which is sufficient for the purpose of collecting statistics.

While spatial averaging along the homogeneous z direction is straightforward, averaging along the streamwise wavy surface can be done using multiple approaches. In this study, we define a local vertical coordinate, d, at each streamwise location with d = 0 at the wall. Its maximum possible value corresponds to the midchannel location and changes with the streamwise coordinate. We then perform streamwise averaging along constant values of d to generate mean statistical profiles. This approach works well for $aδ≪1$ as it tries to approximate the terrain as nearly flat with a large local radius of curvature and therefore nearly homogeneous. To justify this approach, we note that in our study a = 0.07δ (a⁺ ≈ 12, δ = 1), which is an order of magnitude larger than the typical viscous length scale, $Lv$ = ν/u_τ = 1/Re_τ ≈ 0.0055, but smaller than the start of the log layer (y⁺ ≈ 50).

As the mean channel height, δ (for wavy geometry), is kept constant across all the different steepness, ζ, the observed changes in the mean statistics are only due to surface slope effects. In Fig. 6, we show the inner-scaled, double averaged streamwise and vertical velocity for the different surfaces. The different colors in the plot, namely, blue, green, red, lime, orange, and magenta correspond to different wave steepness, ζ = 0, ζ = 0.011, ζ = 0.017, ζ = 0.022, ζ = 0.033, and ζ = 0.044, respectively. The double-averaged, inner-scaled streamwise velocity shows the well-known downward shift of the log-region in the u⁺–y⁺ plot for increasing wave steepness, ζ [Fig. 6(a)], and is indicative of the flow slowing down near the surface from increased drag due to steeper undulations. The double-averaged, inner-scaled vertical velocity structure [Fig. 6(b)] show increasingly stronger net upward flow close to the surface with an increase in ζ. Away from the surface in the logarithmic region, $⟨v⟩x,z,t+$ shows downward flow so that there is no net vertical transport. For the horizontally homogeneous flat channel (ζ = 0), the mean vertical velocity is zero. Therefore, these well-established vertical motions in the mean over wavy surfaces [although small, i.e., ⟨$v$⟩⁺ = O(0.1)] represent the obvious form of surface-induced deviations from equilibrium. Despite these near surface deviations, the dynamics outside the roughness sublayer tend to be similar when normalized and shifted appropriately as illustrated through the defect velocity profiles in Figs. 6(c) and 6(d) that indicate little to no deviation between ζ = 0 and ζ = 0.044.

The normalized mean streamwise velocity gradient helps characterize the different regions of a TBL, especially the inertial (or logarithmic) region. In this study, we estimate the normalized premultiplied inner-scaled mean gradient, $γ=y+d⟨u⟩x,z,t+dy+$ [shown in Fig. 7(a)], which achieves a near constant value of 1/κ (where κ is the von Kármán constant) in the inertial sublayer due to the normalization of the mean gradient by the characteristic law of the wall variables, i.e., surface layer velocity (u_τ) and distance from the wall (y). In this study, for the chosen friction Reynolds numbers Re_τ, we observe that the inertial layer exists over y⁺ ∼ 60–110 for ζ = 0 and shifts to y⁺ ∼ 70–120 for ζ = 0.044, i.e., an upward (rightward) shift in the log layer with the increase in wave steepness, ζ. The estimated von Kármán constants are tabulated in Table II and show a range of 0.39–0.40 for the different runs with κ increasing with ζ. In this study, we use the appropriate value of κ to compute the different metrics. This is in contrast to the flow dependent κ values reported in the work of Leonardi and Castro²⁴ for cube roughness, i.e., κ is reported to decrease from 0.41 in smooth channels to ∼0.35 in the rough-wall TBL.⁵⁸ 

For completeness, we also show the non-dimensional mean streamwise velocity gradient, $Φ=κyuτd⟨u⟩x,z,tdy=γκ$⁠, in Fig. 7(b). The Φ profiles for different ζ mimic the characteristic equilibrium structure starting from zero at the wall followed by a peak at the edge of viscous layer and subsequently a gradual decrease in the buffer layer to a value of one in the inertial sublayer indicative of overall shape similarity in Φ. The shape of the above curves, namely, the upward shift in the log region [Fig. 7(a)] and the smaller peak in Φ with the increase in ζ indicate that the buffer layer becomes increasingly thicker for steeper waves. The “buffer layer” is a region of high turbulence production⁵⁹ where both the viscous and Reynolds stresses are significant. Therefore, the expansion of the buffer layer with ζ is combined with the expansion of the turbulence production zone due to the wavy surface as evidenced from Fig. 8 where the turbulence kinetic energy (TKE) production grows and decays slower for higher ζ in the buffer region (y⁺ ≈ 10–50) in both inner-scaled [Fig. 8(b)] and dimensional [Fig. 8(a)] forms. In fact, this is explicitly seen from the production–dissipation ratio in Fig. 8(c) where the horizontal lines clearly show the upward shift in $⟨Pii⟩x+/⟨Eii⟩x+$ with ζ.

In order to interpret the structure of the components of the Reynolds stress tensor, we need to study its evolution mechanisms. Below, we provide an overview of Reynolds stress transport and the nomenclature adopted over the rest of this manuscript. The Reynolds stress transport equation is shown in Eq. (3). Here, $Lij$ is the local rate of change, and $Cij$ and $Dij$ represent advective and diffusive transport, respectively. The local terms in the evolution equation are $Pij$ representing production, $Eij$ representing dissipation, and $Rij$ is the pressure-rate-of-strain correlation contributing to the redistribution of Reynolds stress. All the above terms are estimated using averaged quantities along the only homogeneous spanwise direction (z) and over a stationary window (t), given by the notation ⟨⟩_z,t. In this study, the stationary window of time is sampled using 2500 temporal snapshots collected over 20 flow through times. Therefore, each of the terms in the following equation varies over the (x, y) space:

In the above, the indices i, j = 1, 2, and 3 correspond to the streamwise (x), vertical (y), and spanwise (z) directions, respectively. In addition, u₁ = u, u₂ = $v$⁠, and u₃ = $w$⁠. On account of statistical stationarity, $Lij=0$⁠, which allows us to rewrite Eq. (3) as

We further average these individual terms along the inhomogeneous streamwise (x) direction, i.e., double averaging. The cumulative effect of the locally acting terms in the transport equation, namely, the sum of production, dissipation, and pressure-rate-of-strain is denoted by $⟨Λij⟩x$⁠, expressed as

Using the above equation, we can indirectly compute the streamwise-averaged diffusive transport term $⟨D⟩x$ as

In Subsections IV B–IV D, we explore the structure of the diagonal elements of the Reynolds stress tensor, $⟨ui′uj′⟩x,z,t,i=j$⁠, with particular focus on the streamwise (double) averaged second order turbulence structure (Fig. 9) and their underlying generation (both double- and single-averaged) in response to the surface undulations.

Noticeable deviations from equilibrium in wavy wall turbulence occur in the second order statistics. We observe in Fig. 9(a) the inner-scaled streamwise variance that peaks in the buffer layer, and this peak value decreases systematically with the increase in ζ. This ζ dependence is localized to the roughness sublayer as the inner scaled profiles collapse in the outer region for all the different cases. The location of peak streamwise variance shifts systematically upward with the increase in ζ, which is enhanced further due to consistent flow separation at larger wave slopes. These observations are consistent with Ref. 35, which shows that the primary influence of effective wave slope, ζ, on the near surface turbulence structure is to accelerate the transition to isotropy of the Reynolds stress tensor as one moves away from the surface. Most literature⁶⁰ suggest the existence of such an upward shift only in response to increases in roughness scale (a⁺) but not so with the effective slope (i.e., λ⁺ for fixed a⁺). This aligns with the classical view of fully rough wall turbulence^17,59 where the peak variances are expected to occur at nearly the roughness height, a. In fact, wall stress boundary conditions for large eddy simulation over rough surfaces are designed to approximate this behavior, especially at high Reynolds numbers with large scale separation (δ/a). We remind the readers to note that, in our study, the wave amplitude, a, is fixed, while the wavelength, λ, is decreased to increase ζ = 2a/λ. For a fixed friction Reynolds number, the decrease in λ (or the increase in ζ) increases the net drag and, in turn, the friction velocity, u_τ. The resulting viscous length scale, $Lv=νuτ$⁠, changes very little (because for this fully developed channel flow, δ and δ⁺ = δ/$Lv$ are constant by design) as does the inner-scaled wave height, a⁺ = a/$Lv$ (which varies from 12.6 to 12.9 for ζ changing from 0.011 to 0.044 as shown in Table I). Therefore, it is evident that this systematic upward shift in the location of peak streamwise variance arises solely from the changes in λ⁺ and ζ. In the following, we explore this in great detail by studying the different terms of the transport equation in Fig. 10.

To dissect the above trends in the inner-scaled double-averaged streamwise variance, we quantify the inner-scaled double-averaged terms in the transport equations (3) and (4), as shown in Fig. 10. For ease of readability, we include the profile of the streamwise variance [Fig. 10(a)] along with profile plots of double-averaged transport equation terms [shown in Figs. 10(b)–10(h)]. Looking at the magnitudes, the dynamics is controlled by the locally operating terms such as production (⁠$⟨P11⟩x+$⁠), dissipation (⁠$⟨E11⟩x+$⁠), and pressure-rate-of-strain (⁠$⟨R11⟩x+$⁠), which are collectively denoted by $⟨Λ11⟩x+$⁠. Of these, we note that production and dissipation are dominant, while the pressure-rate-of-strain term is relatively less important for this streamwise variance transport. Furthermore, the double-averaged cumulative local variance generation rate, $⟨Λ11⟩x+$⁠, matches the double-averaged diffusive transport, $⟨D11⟩x+$⁠. Combining this with Eqs. (4)–(6), we infer that the (double-averaged) advective transport, $⟨C11⟩x+$⁠, has little impact on the evolution of (double-averaged) streamwise variance [see Fig. 23(a) in the  Appendix]. Furthermore, the positive variance production term, $⟨P11⟩x+$⁠, shows the same trend as $⟨u′2⟩x+$ with its peak location trending upwards with the increase in ζ while the peak magnitude decreases. Conversely, the small negative growth rate term, $⟨R11⟩x+$⁠, displays a downward shift in the location of its peak magnitude with ζ so that the effective local turbulence generation, $⟨Λ11⟩x+$⁠, displays an exaggerated upward shift (of its peak). Mechanistically, $⟨R11⟩x+$ represents conversion of streamwise variance to other component variances and is, therefore, negative with peak magnitudes occurring away from the surface due to wall-damping of the vertical turbulent motions. The inner-scaled dissipation rate, $⟨E11⟩x+$⁠, is always positive with magnitudes peaking at the surface and decreasing with ζ through the boundary layer. Taking this order of importance into account among the different terms, we focus primarily on the production structure in the following discussions.

We dissect the inner scaled streamwise variance production [Fig. 10(b)] term, $⟨P11⟩x+$⁠, in the $⟨u′2⟩x,z,t+$ transport equation. In the viscous layer (y⁺ ≲ 7), $⟨P11⟩x+$ is nearly zero and independent of ζ. Beyond this viscous force dominated region, the production increases monotonically to peak in the buffer layer, followed by the subsequent decrease in the log-layer where the different curves collapse (once again independent of ζ). As expected from knowledge of $⟨u′2⟩x,z,t+$ profiles, the magnitude and location of the peak $⟨P11⟩x+$ depends on ζ. We split this component variance production into its dominant contributions, $⟨P11u′u′⟩x+=⟨⟨u′u′⟩z,t+d⟨u⟩z,t+/dx+⟩x$⁠, which can be attributed to surface undulations and $⟨P11u′v′⟩x+=⟨⟨u′v′⟩z,t+d⟨u⟩z,t+/dy+⟩x$ represents production from flow shear, as shown in Figs. 10(c) and 10(d), respectively.

The $⟨u′2⟩x,z,t+$ production arising from the interaction of the inner scaled mean shear, $d⟨u⟩z,t+/dy+$⁠, with the inner scaled vertical momentum flux (⁠$⟨u′v′⟩z,t+$⁠) denoted by $⟨P11u′v′⟩x+$ peaks in the buffer layer [y⁺ ≈ 12–18, as shown in Fig. 10(d)], and these peaks shift systematically upward (see horizontal lines) as ζ varies over 0–0.044. This trend can be interpreted approximately through Figs. 11(a) and 12(b) representing the double averaged profiles of normalized covariance, $⟨u′v′⟩x,z,t+$⁠, and mean shear, $d⟨u⟩x,z,t+/dy+$⁠, respectively. This interpretation is inexact for non-flat wavy surfaces as the streamwise production from the interaction of mean shear with the mean vertical turbulent flux given by $⟨P11u′v′⟩x+=⟨⟨u′v′⟩z,t+d⟨u⟩z,t+/dy+⟩x$ is not the same as $⟨u′v′⟩x,z,t+d⟨u⟩x,z,t+/dy+$⁠. We denote the latter expression as surrogate or pseudo-production, $⟨P11u′v′⟩*+$⁠, which is accurate only in the homogeneity limit.

The double-averaged covariance $⟨u′v′⟩x,z,t+$ peaks at the edge of the buffer layer at y⁺ ≈ 28–34 with the peak height decreasing with ζ [see the horizontal lines in Fig. 11(a)], whereas the normalized mean shear [Fig. 12(b)] achieves its maximum value near the surface. In addition, the peak magnitude of $⟨u′v′⟩x,z,t+$ increases with ζ, while the maximum for $d⟨u⟩x,z,t+/dy+$ decreases. This is consistent with the notion of turbulent stresses forming a higher fraction of the total drag at higher wave slopes, i.e., form drag becomes increasingly important relative to viscous drag, especially for roughness Reynolds numbers, a⁺, being an order of magnitude smaller than the fully rough regime and kept constant (Table I). This aligns with observations for the transitional roughness regime,²² where the surface increasingly moves from the “waviness” to “roughness” regime with an increase in effective slope (2ζ) resulting in higher form drag. The combined influence of the surface-induced trends in $⟨u′v′⟩x,z,t+$ and $d⟨u⟩z,t+/dy+$ [as summarized in Figs. 11(b)–11(d) and 12(b)] show that (i) the Reynolds stress dominates the viscous stresses increasingly closer to the surface at higher ζ and (ii) $⟨P11u′v′⟩x+$ [shown in Fig. 10(d)] shows peak production occurring farther from the surface at higher ζ while decreasing in magnitude. This represents an interesting effect of surface heterogeneity where the peak production does not coincide with the crossover location of viscous and Reynolds stresses as the growth of the latter is steeper compared to a decrease in the former with the height [see Fig. 11(c)].

To decipher the production mechanisms very close to the surface in the viscous layer, we look at Figs. 11(b), 12(b), and 10(d). We see that for ζ sufficiently larger than 0, the vertical turbulent momentum flux $⟨u′v′⟩x,z,t+>0$ [Fig. 11(b)] for y⁺ ≲ 7, resulting in $⟨P11u′v′⟩x+≲0$⁠, i.e., small variance destruction close to the wall. While $⟨u′v′⟩x,z,t+≈0$ close to the surface over flat surfaces (from wall damping), the presence of surface undulations causes larger and positively correlated u′ and $v′$ resulting in ζ dependent destruction of $⟨u′2⟩x,z,t+$⁠. However, the net $⟨u′2⟩x,z,t+$ production, $⟨P11⟩x+$ [in Fig. 10(b)], is nearly zero for y⁺ ≲ 7 with little dependence on ζ due to $⟨P11u′v′⟩x+$ being balanced by surface induced production, $⟨P11u′u′⟩x+=⟨⟨u′2⟩z,t+d⟨u⟩z,t+/dx+⟩x$⁠, i.e., $⟨P11u′u′⟩+>0$ and $⟨P11u′v′⟩+<0$⁠.

We know that $⟨P11u′u′⟩x+=⟨⟨u′2⟩z,t+d⟨u⟩z,t+/dx+⟩x=0$ [shown in Fig. 10(c)] in a flat channel (ζ = 0), whereas for non-flat surfaces (ζ > 0), the streamwise gradient of the mean streamwise velocity (⁠$d⟨u⟩x,z,t+/dx$⁠) is non-zero, resulting in production close to the surface (viscous layer) and its destruction above it in the buffer layer before approaching zero in the logarithmic region. Consequently, there is no significant net turbulence generation over the entire TBL from $⟨P11u′u′⟩x+$⁠, except for pockets of local production and destruction that helps control the shape of overall production $⟨P11⟩x+$⁠. Away from the surface, $⟨P11u′u′⟩x+$ and $⟨P11u′v′⟩x+$ are still complementary, but the different terms do no cancel out as $⟨P11u′v′⟩x+>⟨P11u′u′⟩x+$⁠. Later, we explore whether this is simply a consequence of $d⟨u⟩z,t+/dy+≫d⟨u⟩z,t+/dx+$ (due to ζ ≪ 1) or is more complicated. The different components contribute to the overall production trends as follows: The peak location of $⟨P11u′v′⟩x+$ shows systematic upward shifts with ζ, while its magnitude displays very little sensitivity to the wave slope. In contrast, the profiles of $⟨P11u′u′⟩x+$ show very strong ζ-dependence of the peak destruction magnitudes in the buffer layer, but its location does not. Therefore, $⟨P11⟩x+$ dependence on ζ in both magnitude and shape arises from $⟨P11u′u′⟩x+$ and $⟨P11u′v′⟩x+$⁠, respectively. One can use $⟨P11u′u′⟩x+$ [Fig. 10(c)] to characterize the roughness sublayer height, which, in this case, is ∼3a⁺ and nearly independent of ζ.

$P11+(x,y)$ represents the inner-scaled variance production based on averaging along the homogeneous direction (z) and over a stationary window (t) of the turbulent flow. $P11+(x,y)$ shown in Fig. 13 is averaged along the inhomogeneous streamwise (x) direction to estimate the averaged production, $⟨P11⟩x+$⁠, shown in Fig. 10(b). Interpretation of $⟨P11⟩x+$ requires understanding the structure of $P11+(x,y)$ and its components. Figure 13 shows the inner-scaled production over the y–ϕ space, with ϕ being the non-dimensional streamwise phase coordinate (ϕ = 2πx/λ) and y = y/δ for unit half channel height (δ = 1). We clearly see that the structure of averaged production contours is qualitatively invariant in this y–ϕ space for different ζ = 2a/λ > 0 while the magnitude depends on the wave slope with higher ζ producing stronger peaks and troughs. To identify the different layers of the TBL, we define d as the local vertical coordinate relative to the wall at each streamwise location.

Looking at the isocontours in Figs. 13(b) and 13(c), we note that both $P11u′v′+(x,y)$ and $P11u′u′+(x,y)$ indeed play complementary roles of production and destruction in different regions of TBL, especially close to the surface (d ≲ 0.1). Specifically, $P11u′v′+(x,y)<0$ (destruction) along the windward side of the wavy surface (i.e., ϕ = 0 − π and d ≲ 0.1) as indicated by the cyan region in Fig. 13(c), while along the leeward side (i.e., ϕ = π − 2π and d ≲ 0.1), there exists very little turbulence production (shown in black). The exception to this is a small negative production (cyan) zone in the trough due to flow separation at larger ζ. This is complemented by $P11u′u′+(x,y)>0$ [in Fig. 13(b)] along the windward side (i.e., ϕ ≈ 0 − π and d ≲ 0.1) and near-zero values along the leeward side (i.e., ϕ ≈ π − 2π and d ≲ 0.1). This complementary structure of $P11u′u′+$ and $P11u′v′+$ when averaged along the streamwise direction yield net production, $⟨P11u′u′⟩x+>0$ [Fig. 10(c)], and destruction, $⟨P11u′v′⟩x+<0$ [Fig. 10(d)], respectively. Just as $⟨P11u′u′⟩x+$ and $⟨P11u′v′⟩x+$ cancel each other in the inner layer, $P11u′v′+(x,y)$ and $P11u′u′+(x,y)$ also show overlapping regions of production and destruction so that there is no net variance generation in the viscous dominated region over the surface for all ζ [see $P11+(x,y)$ in Fig. 13(a)]. Away from the surface in the outer layer, $P11u′v′+(x,y)$ [Fig. 13(c)] dominates and controls the large-scale structure of $P11+(x,y)$⁠. To interpret this structure, the following observations were made regarding the relevant strain rate and Reynolds stress tensor terms:

• On average, $‖d⟨u⟩z,t+/dx+‖≪‖d⟨u⟩z,t+/dy+‖$ by a factor of ∼O(20) [as shown in Figs. 14(b) and 14(d)], which is comparable to O(1/ζ).

• $‖⟨u′2⟩z,t+‖>‖⟨u′v′⟩z,t+‖$ by a factor of ∼O(5) [in Figs. 14(a)–14(c)] and achieve their maximum magnitudes along the surface in the buffer layer.

• It is well known that as streamlines wrap around the wave crest, the flow accelerates creating a local low pressure region over the hump. This shape effect on the streamwise velocity field is felt away from the surface.

• $d⟨u⟩z,t+/dx+$ has an asymmetric structure [Fig. 14(b)] that is different away and close to the surface. Away from the surface, the accelerating and decelerating flow before and after the crest results in $d⟨u⟩z,t+/dx+>0$ and $d⟨u⟩z,t+/dx+<0$⁠, respectively. In the viscous region, this trend is reversed due to the surface slope-induced blockage causing $d⟨u⟩z,t+/dx+<0$ and $d⟨u⟩z,t+/dx+>0$ in the windward and leeward sides.

• $d⟨u⟩z,t+/dy+$ originates primarily from flow shear and, therefore, is positive all along the surface while decreasing rapidly with the height [Fig. 14(d)]. The exception being flow separation at large enough ζ that causes $d⟨u⟩z,t+/dy+<0$ near the trough.

• While the magnitudes of $d⟨u⟩z,t+/dx+$ and $d⟨u⟩z,t+/dy+$ are large over a thin layer along the windward side, we see a more diffused layer {thickness [O(a⁺)]} along the leeward side due to wake mixing behind the crest.

Obviously, variance production/destruction is large in regions of strong correlation between the appropriate component of Reynolds stress and strain rate tensors. Given that $d⟨u⟩z,t+/dx+$⁠, $d⟨u⟩z,t+/dy+$ [Figs. 14(b)–14(d)] are strong near the surface while $⟨u′2⟩z,t+$⁠, $⟨u′v′⟩z,t+$ achieve larger magnitudes away from the wall [Figs. 14(a)–14(c)], the strong destruction/production zone in both $P11u′u′+$ and $P11u′v′+$ exists [Figs. 13(a)–13(b)] in the lower buffer region. However, as $d⟨u⟩z,t+/dx+$ and $d⟨u⟩z,t+/dy+$ assume a more diffused structure behind the wave crest due to wake mixing, the resulting production/destruction zone is also thicker with larger magnitudes in this leeward region and grows with slope, ζ. For such wavy surfaces, the leeward side production (as well for the windward side away from the surface) is negative for $P11u′u′+$ and positive for $P11u′v′+$⁠. Given that $‖P11u′u′+‖<‖P11u′v′+‖$ over most of the TBL, we see that the structure of net production, $P11+$⁠, is dominated by $P11u′v′+$ that depends on $d⟨u⟩z,t+/dy+$ [Fig. 14(d)] and $⟨u′v′⟩z,t+$ [Fig. 14(c)]. Specifically, the primary generation of $⟨u′2⟩z,t+$ [see $P11+$ in Fig. 13(a)] occurs in the thick wake-induced buffer region along the leeward slope through shear production, $P11u′v′+$⁠. However, the windward slope is responsible for surface induced near surface variance generation through $P11u′u′+$ [i.e., $d⟨u⟩z,t+/dx+$ Fig. 14(b)] to overcome the destruction contained in $P11u′v′+$ due to positive covariance, $⟨u′v′⟩z,t+>0$ [cyan region in Fig. 14(c)]. In addition to increasing strain rates at higher ζ, the flow also involves unsteady separation dynamics with $d⟨u⟩z,t+/dy+<0$ [see cyan regions near the wave trough shown in Fig. 14(d)] and turbulence destruction zones [see cyan region in Figs. 13(a) and 13(c) for ζ = 0.044] that grow thicker with ζ. These results suggest that the influence of ζ shows up in at least three ways: (i) larger mixing and mean shear in the leeward side of the wave resulting in enhanced production, (ii) surface induced near-wall variance generation along the windward side, and (iii) separation-induced destruction in the trough.

The two-dimensional surface undulations generate a complex production structure in $P11+(x,y)$ (Fig. 13) for ζ > 0 that is submerged within its one-dimensional surrogate, $⟨P11⟩x+$ [Fig. 10(b)]. There exist multiple ways to characterize the surface dispersion effects on turbulent statistics⁶¹ in order to estimate the roughness sublayer height. Here, we characterize the surface dispersion effects using the surrogate or pseudo-production estimate, which computes the Reynolds stress–strain rate interaction as if homogeneity exists, i.e., as the product of the double averaged Reynolds stress tensor (i.e., $⟨u′u′⟩x,z,t+,⟨u′v′⟩z,t+$⁠) and the mean strain rate tensor (⁠$d⟨u⟩x,z,t+/dy+,d⟨u⟩x,z,t+/dx+,…$⁠). The formal error contained in this production estimate given by $⟨P11u′v′⟩x+−⟨P11u′v′⟩*+=⟨⟨u′v′⟩z,t+d⟨u⟩z,t+/dy+⟩x−⟨u′v′⟩x,z,t+d⟨u⟩x,z,t+/dy+$ represents surface dispersion effects contained in the computation of turbulence production. This is quantified in Fig. 15 with the top row representing the pseudo-estimates and the bottom row representing the deviations. Obviously, in the homogeneous limit (ζ = 0), the dispersion in production is zero while it grows systematically with wave slope, ζ. The dispersion errors are concentrated closer to the surface (i.e., y⁺ ≲ 50 ∼ 4a⁺) and decrease gradually to zero in the outer region. The surface influence on the TBL (∼4a⁺) extends further than that observed for the surface-induced production, $⟨P11u′u′⟩x+$ (∼3a⁺). While the different pseudo-production estimates underpredict and overpredict depending on the region of the TBL, it shares a qualitative similarity with the true estimates.