Direct numerical simulation of turbulence over two-dimensional waves

Knowledge of turbulent flows over non-flat surfaces is of major practical interest in diverse applications. Significant work continues to be reported in the roughness regime at high Reynolds numbers where the cumulative effect of surface undulations on the averaged and integrated turbulence quantities is well documented. Even for such cases, the surface topology plays an important role for transitional roughness Reynolds numbers that is hard to characterize. In this work, we explore in detail the mechanisms underlying turbulence generation and transport, particularly within the region of the turbulent boundary layer affected by the surface. We relate surface shape with turbulence generation mechanisms and Reynolds stress transport, which has implications to drag increase. We accomplish this using a suite of direct numerical simulations of fully developed turbulent flow between two infinitely wide, two-dimensional sinusoidally wavy surfaces at a friction Reynolds number, Reτ = 180, with different mean surface slopes, ζ (and fixed inner-scaled undulation height, a+ = 13), corresponding to the “waviness” regime. The increase in wave slope enhances near surface turbulent mixing, resulting in increased total drag, higher fraction of form drag, faster approach to isotropy, and thereby modulation of the buffer layer. The primary near-surface streamwise and vertical turbulence production occurs in the leeward and windward sides of the wave, respectively. Furthermore, this production structure shows significant dispersion effects. However, the primary generation process of vertical and spanwise variance occurs through pressure–strain rate mechanism in the windward side.Knowledge of turbulent flows over non-flat surfaces is of major practical interest in diverse applications. Significant work continues to be reported in the roughness regime at high Reynolds numbers where the cumulative effect of surface undulations on the averaged and integrated turbulence quantities is well documented. Even for such cases, the surface topology plays an important role for transitional roughness Reynolds numbers that is hard to characterize. In this work, we explore in detail the mechanisms underlying turbulence generation and transport, particularly within the region of the turbulent boundary layer affected by the surface. We relate surface shape with turbulence generation mechanisms and Reynolds stress transport, which has implications to drag increase. We accomplish this using a suite of direct numerical simulations of fully developed turbulent flow between two infinitely wide, two-dimensional sinusoidally wavy surfaces at a friction Reynolds number, Reτ = 180, with different mean surf...


I. INTRODUCTION
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, 1 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 2 nearly two hundred years ago; in the early half of the 20th century by Nikuradse, 3 Colebrook, 4 and Moody; 5 and more recently by various research groups. [6][7][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][10][11][12][13][14][15][16] and reviewed in Refs. 14 and 17. In addition, there exist studies of flow over systematically designed roughness using experiments 9 and more recently using direct numerical simulation (DNS) 8,[18][19][20][21][22][23][24][25][26][27][28] and large eddy simulation (LES). 29 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. 32 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][34][35] and elsewhere. 36 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.

A. Relevance to roughness characterization
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, 39 Δ⟨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.

B. Turbulence structure over wavy surface undulations
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][43][44][45][46][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 48 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 49 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 18 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 20 showed that surface undulations can lead to turbulence generation even for linearized flow dynamics. Wavy surface undulations impact the nearwall 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 50 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 secondorder turbulence structure, namely, the components of the Reynolds stress tensor and their generation mechanisms as modulated by the

II. NUMERICAL METHODS
In this study, we adopt a customized version of the Incom-pact3D 50 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. 53 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 52 differ in the details of this reconstruction. In this study, we adopt the one-dimensional higher order polynomial reconstruction of Gautier, Laizet, and Lamballais 54 and refer to the work of Khan and Jayaraman 35 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 26 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.

A. Simulation design
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, Re b = u b δ ν , 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 55 and comparable to the full span design of MacDonald et al. 25 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. 26 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 twodimensional undulations. Given the higher numerical accuracy of our method compared to the DNS algorithms of MacDonald et al. 26 and Jelly and Busse 55 who used 2nd order schemes as compared to TABLE I. Tabulation of different design parameters for the simulations such as wavelength (λ), amplitude (a) and steepness (ζ = 2a λ ) of the wavy surface, friction velocity (uτ), Reynolds numbers (Re) based on the boundary layer height (δ), and different velocities expressed as subscripts ("cl" = centerline quantity, "b" = bulk quantity, and "τ" = based on friction velocity) and the grid spacing in different directions ("Δx" = streamwise grid spacing, "Δz" = spanwise grid spacing, "Δyw" = wall normal grid spacing near the surface, "Δy cl " = wall normal grid spacing near the flow centerline, and "Δt" = simulation time step). The superscript "+" refers to the inner scaled quantity [scaled with respect to kinematic viscosity (ν) and friction velocity (uτ)]. 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.

B. Convergence of turbulence statistics
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), 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 ⟩ + , 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.

C. Assessment of simulation accuracy
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 56 regenerated by Moser, Kim, and Mansour 57 (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 6 (128 × 129 × 128) grid points and solved the flow equations by advancing modified variables, namely, wallnormal 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

III. MEAN FIRST-ORDER TURBULENCE STRUCTURE
The primary goal of this study is to explore the nonequilibrium, 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 (streamwiseaveraged) 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 streamwiseaveraged 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)

A. Streamwise averaging of turbulence 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).

B. Outer layer similarity and mean velocity profiles
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, ζ 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.

C. Quantification of mean velocity gradients and inertial sublayer
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⟩ +  Table I).

FIG. 7.
Variation of non-dimensional mean streamwise velocity gradients: (a) 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 24 for cube roughness, i.e., κ is reported to decrease from 0.41 in smooth channels to ∼0.35 in the rough-wall TBL. 58 For completeness, we also show the non-dimensional mean streamwise velocity gradient, Φ = κy u τ d⟨u⟩ x,z,t dy = γκ, 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 59 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 ζ.

A. Overview of Reynolds stress transport
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 1 = u, u 2 = v, and u 3 = 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 60 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 , and diffusion (h). We further split the production term ⟨P11⟩ + The horizontal lines correspond to the vertical location of the maximum/minimum value for a chosen statistic. If the peak locations are different, we color match the horizontal lines with the corresponding curves.

ARTICLE
scitation.org/journal/adv 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.

Dynamics of streamwise variance transport
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 Looking at the magnitudes, the dynamics is controlled by the locally operating terms such as production (⟨P 11 ⟩ + x ), dissipation (⟨E 11 ⟩ + x ), and pressurerate-of-strain (⟨R 11 ⟩ + 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 doubleaveraged cumulative local variance generation rate, ⟨Λ 11 ⟩ + x , matches the double-averaged diffusive transport, ⟨D 11 ⟩ + x . Combining this with Eqs. (4)-(6), we infer that the (double-averaged) advective transport, ⟨C 11 ⟩ + x , has little impact on the evolution of (doubleaveraged) streamwise variance [see Fig. 23(a) in the Appendix]. Furthermore, the positive variance production term, ⟨P 11 ⟩ + 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, ⟨R 11 ⟩ + 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, ⟨R 11 ⟩ + 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 innerscaled dissipation rate, ⟨E 11 ⟩ + 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.

Mechanisms of streamwise variance generation
We dissect the inner scaled streamwise variance production [ Fig. 10(b)] term, ⟨P 11 ⟩ + x , in the ⟨u ′2 ⟩ + x,z,t transport equation. In the viscous layer (y + ≲ 7), ⟨P 11 ⟩ + 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 ⟨P 11 ⟩ + x depends on ζ. We split this component variance production into its dominant contributions, ⟩ + 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 ( 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 x,z,t /dy + . We denote the latter expression as surrogate or pseudo-production, ⟨P u ′ v ′ 11 ⟩ + * , 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, 22 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) 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 ⟨P u ′ v ′ 11 ⟩ + 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, ⟨P 11 ⟩ + x [in Fig. 10(b)], is nearly zero for y + ≲ 7 with little dependence on ζ due to ⟨P u ′ v ′ 11 ⟩ + x being balanced by surface induced production, 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 ⟨P u ′ u ′ 11 ⟩ + x , except for pockets of local production and destruction that helps control the shape of overall production ⟨P 11 ⟩ + x . Away from the surface, ⟨P u ′ u ′ 11 ⟩ + x and ⟨P u ′ v ′ 11 ⟩ + x are still complementary, but the different terms do no cancel out as ⟨P u ′ v ′ 11 ⟩ + x > ⟨P u ′ u ′ 11 ⟩ + 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 ⟨P u ′ v ′ 11 ⟩ + x shows systematic upward shifts with ζ, while its magnitude displays very little sensitivity to the wave slope. In contrast, the profiles of ⟨P u ′ u ′ 11 ⟩ + x show very strong ζ-dependence of the peak destruction magnitudes in the buffer layer, but its location does not. Therefore, ⟨P 11 ⟩ + x dependence on ζ in both magnitude and shape arises from Fig. 10(c)] to characterize the roughness sublayer height, which, in this case, is ∼3a + and nearly independent of ζ. c. Two-dimensional structure of ⟨u ′2 ⟩ + z,t production. P + 11 (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. P + 11 (x, y) shown in Fig. 13 is averaged along the inhomogeneous streamwise (x) direction to estimate the averaged production, ⟨P 11 ⟩ + x , shown in Fig. 10(b). Interpretation of ⟨P 11 ⟩ + x requires understanding the structure of P + 11 (x, y) and ARTICLE scitation.org/journal/adv FIG. 13. Contours of the inner-scaled (single) averaged streamwise variance production, P + 11 (a) and its components P u ′ u ′ 11 + (b) and P u ′ v ′ 11 + (c) as a function of ϕ and y, where ϕ = 2πx/λ. The cyan colored region represents negative production of streamwise variance.
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 P u ′ v ′ 11 + (x, y) and P u ′ u ′ 11 + (x, y) indeed play complementary roles of production and destruction in different regions of TBL, especially close to the surface (d ≲ 0.1). Specifically, P u ′ v ′ 11 + (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 P u ′ u ′ 11 + (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 P u ′ u ′ 11 + and P u ′ v ′ 11 + when averaged along the streamwise direction yield net production, ⟨P u ′ u ′ 11 ⟩ + x > 0 [ Fig. 10(c)], and destruction, ⟨P u ′ v ′ 11 ⟩ + x < 0 [ Fig. 10(d)], respectively. Just as ⟨P u ′ u ′ 11 ⟩ + x and ⟨P u ′ v ′ 11 ⟩ + x cancel each other in the inner layer, P u ′ v ′ 11 + (x, y) and P u ′ u ′ 11 + (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 P 11 + (x, y) in Fig. 13(a)]. Away from the surface in the outer layer, P u ′ v ′ 11 + (x, y) [ Fig. 13(c)] dominates and controls the largescale structure of P 11 + (x, y). To interpret this structure, the following observations were made regarding the relevant strain rate and Reynolds stress tensor terms: i. 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 compa-  iii. 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. iv. 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. v. 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. vi. 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 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 P u ′ u ′ 11 + and positive for P u ′ v ′ 11 + . Given that ∥P u ′ u ′ 11 + ∥ < ∥P u ′ v ′ 11 + ∥ over most of the TBL, we see that the structure of net production, P 11 + , is dominated by P u ′ v ′ 11 + that depends on d⟨u⟩ + z,t /dy + [ Fig. 14(d)] and Fig. 14(c)]. Specifically, the primary generation of ⟨u ′2 ⟩ + z,t [see P 11 + in Fig. 13(a)] occurs in the thick wake-induced buffer region along the leeward slope through shear production, P u ′ v ′ 11 + . However, the windward slope is responsible for surface induced near surface variance generation through P u ′ u ′ 11 + [i.e., d⟨u⟩ + z,t /dx + Fig. 14(b)] to overcome the destruction contained in P u ′ v ′ 11 + 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 ARTICLE scitation.org/journal/adv 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.
d. Dispersion effects in production. The two-dimensional surface undulations generate a complex production structure in P + 11 (x, y) (Fig. 13) for ζ > 0 that is submerged within its one-dimensional surrogate, ⟨P 11 ⟩ + x [ Fig. 10(b)]. There exist multiple ways to characterize the surface dispersion effects on turbulent statistics 61 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 stressstrain 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 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

ARTICLE scitation.org/journal/adv
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, ⟨P u ′ u ′ 11 ⟩ + 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.
The effect of surface undulations with increasing slope, ζ, on vertical variance profiles [ Fig. 16(a)] is to enhance its growth near the surface, especially in the buffer layer, i.e., y + ≈ 7-40. The steeper surface undulations increase the fraction of the vertical variance contributing to the turbulent kinetic energy (TKE) as was observed in Ref. 35. It is well known that the vertical variance near the surface is , and diffusion (h). We further split the production term ⟨P22⟩ + The horizontal lines correspond to the vertical location of the maximum/minimum value for a chosen statistic. If the peak locations are different, we color match the horizontal lines with the corresponding curves.

ARTICLE
scitation.org/journal/adv smaller compared to the streamwise variance due to the wall damping effect, which, in turn, causes the variance to peak further away from the surface. For these cases, the peak vertical variance occurs in the buffer-log transition region (y + ≈ 50-55 ≳ 4a + ), which is outside the roughness sublayer (i.e., ≈3a + ). Therefore, it is not surprising that the peak vertical variance magnitude and location is relatively insensitive to ζ as farther away from the surface (beyond the roughness sublayer), the wall effects have died down. In this region, the turbulence structure is insensitive to ζ with all the different curves collapsing on top of each other and thereby supporting the outer layer similarity.

Dynamics of vertical variance transport
In a TBL over a flat surface, the production of vertical variance from the interaction of the mean strain rate with the Reynolds stresses, ⟩ + x,z,t is generated through a conversion process of streamwise turbulent fluctuations through the pressure-rate-of-strain term, ⟩ + x,z,t is controlled by surface induced production along with mechanisms such as return-to-isotropy, dissipation, and diffusion. The relative magnitudes of the different terms in the variance transport equation show that production (⟨P 22 ⟩ + x ), dissipation (⟨E 22 ⟩ + x ), and pressurerate-of-strain (⟨R 22 ⟩ + x ) as cumulatively depicted by ⟨Λ 22 ⟩ + x dominate variance generation with the production being the smallest. This local generation of vertical variance is balanced only by diffusive transport of the double-averaged, inner-scaled variance, ⟨D 22 ⟩ + x , as the averaged advective transport is insignificant [see Fig. 23(b) in the Appendix] for this statistically stationary system. The larger vertical variance in the buffer layer for higher ζ [ Fig. 16(a)] is also observed in production, ⟨P 22 ⟩ + x [ Fig. 16 In the following, we delve into the mechanisms underlying vertical variance generation including the relatively small, but dynamically important, surface-induced production.

Mechanisms of vertical variance generation
For ζ > 0, the fundamental modulations to ⟨v ′2 ⟩ + x,z,t production occur outside the viscous region of the roughness sublayer as seen from the inner-scaled, averaged vertical variance production, Fig. 16(b). The extent of this production deviates farther from zero with the increase in wave steepness, ζ, as both the streamwise (⟨d⟨v⟩z,t/dx⟩ + x ) and vertical (⟨d⟨v⟩z,t/dy⟩ + x ) gradients of the mean vertical velocity [see Figs. 12(c) and 12(d)] increasingly depart from zero due to local heterogeneity. Consequently, turbulence production, ⟨P v ′ u ′ 22 ⟩ + x = −⟨⟨v ′ u ′ ⟩z,td⟨v⟩z,t/dx⟩ + x , generates vertical variance in the buffer layer and destroys some of it in the viscous layer below [Fig. 16(c)]. Similarly, destroys variance in the viscous layer and generates above it [ Fig. 16(d)]. The ⟨v ′2 ⟩ + x,z,t production, Figs. 16(b) and 16(c)]. This is consistent with the trends in magnitudes of ⟨d⟨v⟩z,t/dx⟩ + x and ⟨d⟨v⟩z,t/dy⟩ + x [see Figs. 12(c) and 12(d)], i.e., ⟨d⟨v⟩z,t/dy⟩ + This suggests that the production process depends strongly on gradients of vertical velocity although the precise nature of this relationship needs to be explored.
In addition to the surface-induced production, the dominant contribution to the local generation of ⟨v ′2 ⟩ + x,z,t , ⟨Λ 22 ⟩ + x comes from ⟨R 22 ⟩ + x . For these small values of ζ, ⟨R 22 ⟩ + x and ⟨λ 22 ⟩ + x are similar in structure with those for homogeneous flat channel turbulence. We observe that closer to the surface in the viscous layer, ⟨R 22 ⟩ + x < 0 due to splat events from wall blockage, i.e., ⟨v ′2 ⟩ + x,z,t is converted to ⟨u ′2 ⟩ + x,z,t and ⟨w ′2 ⟩ + x,z,t . Away from the surface in the buffer layer, x > 0 enables return to isotropy. The increase in ζ enhances this pressure-rate-of-strain mechanism resulting in faster growth of the vertical variance ⟨v ′2 ⟩ + x,z,t through the buffer layer [ Fig. 9(b)] and return to isotropy as evidenced by the location of peak ⟨R 22 ⟩ + x moving closer to the surface.
a. Averaged production of surface undulations. ⟨P 22 ⟩ + x originates solely from heterogeneity effects through Figs. 16(c) and 16(d), respectively. Unlike ⟨P 11 ⟩ + x , both these components for ⟨P 22 ⟩ + x produce vertical variance in the buffer layer (i.e., reinforce each other) due to non-zero streamwise and vertical gradients of ⟨v⟩ + z,t [see Figs. 12(c) and 12(d)]. In fact, ζ directly enters ⟨P 22 ⟩ + x as it relates to the ratio of ⟨d⟨v⟩z,t/dx⟩ + x and ⟨d⟨v⟩z,t/dy⟩ + x . Given that (i) ⟨v ′ u ′ ⟩x,z,t < 0 [everywhere except for a small region near the surface at higher ζ as shown in Fig. 11(b)] and ⟨v ′ v ′ ⟩x,z,t > 0 across the TBL and (ii) ⟨d⟨v⟩z,t/dx⟩ + x ≳ 0 and ⟨d⟨v⟩z,t/dy⟩ + x < 0 in the buffer layer, one can approximately estimate the production These deviations represent dispersion effects on production given by , which are both zero in the absence of surface heterogeneity. Figure 17 shows the pseudoproduction estimates in the top row and the corresponding dispersion contribution in the bottom row. The extent of deviations only increase with ζ. Naturally, both the surface-induced production can be used to characterize the roughness sublayer height. In this case, the average estimate is ∼4a + ≈ 50 for the former and ∼5a + ≈ 60 in the latter with no visible dependence on ζ.
b. Two-dimensional structure of ⟨v ′2 ⟩ + z,t production. For a deeper understanding of the production mechanisms, we look at the inner-scaled two-dimensional production contours of P + 22 (x, y), + (x, y) based on averaging along the homogeneous direction (z) and over a stationary window (t) of the turbulent flow in Fig. 18. We clearly see from the isocontours [Figs. 18(a) and 18(b)] that qualitatively P 22 , and this primary contribution originates from the interaction of vertical variance, ⟨v ′ v ′ ⟩ + z,t , with the vertical gradient of ⟨v ′ ⟩ + z,t . The secondary contribution, P v ′ u ′ 22 + , arises from the interaction of vertical turbulent momentum flux, ⟨v ′ u ′ ⟩ + z,t , with the streamwise gradient of ⟨v ′ ⟩ + z,t . Figure 18 shows the inner-scaled production over the y − ϕ space, where ϕ is a non-dimensional streamwise phase coordinate (ϕ = 2πx/λ) and y = y/δ (δ = 1). We clearly observe that the production contour shape is mostly qualitatively invariant with respect to the wave shape for different ζ = 2a/λ > 0, indicating that the size of the structures scale with λ. However, the inner-scaled magnitude of production (i.e., the colors in Fig. 18) depends on the wave slope with higher ζ producing stronger peaks and troughs. As before, we characterize the regions near and far from the surface using the wall-normal distance relative to the local surface height, d.
Looking at the structure of P + 22 (x, y), we see that the bulk of ⟨v ′2 ⟩ + z,t production from surface undulations occur along the slopes of the wave with a slight departure from the wall (d/δ ≳ 0.05) due to wall blockage. Visually, this structure of P + 22 is derived primarily from P v ′ v ′ 22 + with a secondary contribution from P v ′ u ′ 22 + . Correlating Fig. 18 is characterized by d⟨v⟩ + z,t /dy + . Given this dependence on the strain rate structure, it is worth looking at the structure of mean vertical velocity, ⟨v⟩z"t, generated by the streamlines curving over the wavy undulations. This creates a region of upward flow (⟨v⟩z,t > 0) over the windward slope and downward flow (⟨v⟩z,t < 0) along the leeward slope. These vertical flow structures mostly represent surface form/shape influences and, therefore, extend through most of the roughness layer. This is also true for d⟨v⟩ + z,t /dx + [see Fig. 19(b)], which varies gradually along the vertical direction. The larger gradient, d⟨v⟩ + z,t /dy + (due to small ζ), mostly represents near-wall shear stress and, therefore, decreases rapidly away from the surface [see Fig. 19(d)]. Therefore, d⟨v⟩ + z,t /dx + , despite being smaller in magnitude, persists farther away from the surface. Although it is in the buffer and inertial regions that occurs close to the surface along the leeward slope and away from the surface in the windward side. The smaller production from P u ′ v ′ 22 + occurs away from the surface along the windward slope of the wave to supplement that from P v ′ v ′ 22 + . In summary, the production of ⟨v ′2 ⟩ + z,t dominates in the buffer layer over the windward slope and closer to the surface along the leeward slope. Given the preponderance on strain rates, the higher slopes naturally generate stronger variance production.
D. Spanwise variance, ⟨w ′2 ⟩ + x,z,t Similar to ⟨v ′2 ⟩ + x,z,t , the inner-scaled spanwise variance, ⟨w ′2 ⟩ + x,z,t , also shows stronger growth [see Figs. 9(c) and 20(a)] through the viscous and lower regions of the buffer layer with steeper surface undulations. Additionally, ⟨w ′2 ⟩ + x,z,t peaks in the buffer layer (y + ≈ 35), and this peak location is lower than that observed for ⟨v ′2 ⟩ + x,z,t due to faster growth near the surface (see Fig. 9) in the absence of wall blockage. Furthermore, the peak magnitudes of ⟨w ′2 ⟩ + x,z,t increase with ζ, while the locations of the peak decrease with ζ suggesting surface-enhanced return to isotropy. Unlike the trends for ⟨v ′2 ⟩ + x,z,t , the peak magnitudes and locations for ⟨w ′2 ⟩ + x,z,t show increased sensitivity to ζ as they occur well within the roughness sublayer. Particularly, we see two clusters of peak locations corresponding to ζ = 0.0-0.022 and ζ = 0.033-0.044. This suggests possible dependence of ⟨w ′2 ⟩ + x,z,t on case specific flow features such as flow separation at higher ζ, which, in turn, depends on the Reynolds number. However, the variation in peak magnitudes of ⟨w ′2 ⟩ + x,z,t with ζ is systematic.
Reynolds stress-strain rate interactions [⟨P 33 ⟩ + x = 0 in Fig. 20(b)] due to spanwise homogeneity. Nevertheless, the transport process is still sensitive to ζ due to coupling across the different velocity components. For example, wall-blockage converts vertical motions into other components for all ζ. In fact, we clearly see ⟨R 11 ⟩ + x ≲ 0, ⟨R 22 ⟩ + x ≲ 0, and ⟨R 33 ⟩ + x > 0 for ζ ≥ 0 near the surface in Figs. 10(f), 16(f), and 20(d), respectively. Away from the surface, we have ⟨R 11 ⟩ + x < 0, ⟨R 22 ⟩ + x > 0, and ⟨R 33 ⟩ + x > 0. This dynamics is governed by ⟨R 11 ⟩ + x + ⟨R 22 ⟩ + x + ⟨R 33 ⟩ + x = 0 for incompressible flow, which leaves ⟨R 33 ⟩ + x > 0 (and ⟨R 11 ⟩ + x ≤ 0) through the TBL. In this study, ⟨w ′2 ⟩ + x,z,t is generated closer to the surface using two different pressure-rate-of-strain mechanisms: (i) conversion of diffused ⟨v ′2 ⟩ + x,z,t and (ii) conversion of ⟨u ′2 ⟩ + x,z,t produced by the interaction of mean shear with Reynolds stress. With the increase in effective wave-slope, ζ, ⟨R 11 ⟩ + x , ⟨R 22 ⟩ + x , and ⟨R 33 ⟩ + x not only increase in magnitude through the roughness layer but also assume peak values closer to the surface indicating return to isotropy nearer to the wall. The effect of surface undulations on the different ⟨R⟩ + x profiles is felt until y + = 3a + .
On the windward side before the crest (ϕ ≈ π/4 − π for ζ > 0), R + 11 < 0 (and this magnitude increases with ζ), indicating that splattype events are less likely in this region. This region corresponds , and diffusion (f). The horizontal lines correspond to the vertical location of the maximum/minimum value for a chosen statistic.
to a favorable pressure gradient zone with ⟨u ′ v ′ ⟩ + z,t > 0 [ Fig. 19(a)] associated with surface-induced ⟨u ′2 ⟩ + z,t generation. The inner-scaled pressure-rate-of-strain terms, here, are such that R + 11 < 0, R + 22 < 0, and R + 33 > 0, indicating conversion of ⟨u ′2 ⟩ + z,t and ⟨v ′2 ⟩ + z,t (from diffusion) to ⟨w ′2 ⟩ + z,t . Naturally, this is a region of strong spanwise variance generation [see Figs. 21(a)-21(c)], which is enhanced further by ζ. At sufficiently large ζ, the flow behind the crest (π < ϕ ≲ 3π/2) is impacted by the surface curvature through the pressure field in a way, which pushes the splat dynamics further downstream (closer to the trough) causing R + 11 to be weakly negative. Instead, in this region, the pressure-rate-of-strain terms tend to be small with R + 11 ≲ 0, R + 33 ≳ 0, and R + 22 ≲ 0 or ≳ 0 depending on ζ. This is not surprising given the very little ⟨v ′2 ⟩ + z,t generation (see P + 22 ≈ 0 in Fig. 18) as compared to significant ⟨u ′2 ⟩ + z,t generation (P + 11 > 0 in Fig. 13) over this region. Thus, we end up with three distinct regions along the wavy surface with different pressure-rate-of-strain dynamics.
Away from the surface, the pressure-rate-of-strain term pushes the flow toward isotropy with the dominant ⟨u ′2 ⟩ + z,t increasingly converted to ⟨v ′2 ⟩ + z,t and ⟨w ′2 ⟩ + z,t such that R + 11 < 0, R + 22 > 0, and R + 33 > 0. In fact, the regions with large magnitudes of R + 11 and R +

FIG. 22.
Vertical variation of the ratio of different diagonal elements of averaged pressure-rate-of-strain tensor corresponding to streamwise (a), vertical (b), and spanwise (c) variance. increase. This suggests more rapid pressurerate-of-strain dynamics despite the smaller normalized variances at higher ζ. In fact, this return to isotropy occurs much closer to the mean surface height at higher ζ as evidenced by the faster approach of the ratio R + 33 /R + 22 to unity [ Fig. 22(a)] at higher ζ. These ratios of double-averaged pressure-rate-of-strain terms show that sufficiently far away from the surface, rates of conversion of energy of u ′ fluctuations to that of v ′ and w ′ are equal, i.e., R + 33 /R + 22 ∼ 1 [ Fig. 22(a)

V. CONCLUSION
We present a DNS-based study of turbulence structure over non-flat surfaces, with emphasis on diagonal components of the Reynolds stress and terms that govern their evolution, especially in the region of the TBL affected by surface heterogeneity. For this reason, the high-fidelity DNS is carried out at smaller turbulent Reynolds numbers (Reτ = 180) between two infinitely parallel plates with the wavy roughness mirroring each other. We characterize the shape of different wavy surfaces using an effective slope measure denoted by ζ = 2a/λ. Fixing the roughness Reynolds number and wave amplitude a, we vary ζ over 0-0.044 corresponding to mildly separated flows. Consistent with the literature on rough wall turbulence, the streamwise mean velocity structure indicates a characteristic downward shift of the logarithmic region indicating increased flow drag with wave slope, ζ. This is accompanied by sustained upward vertical flow in the lower roughness sublayer and downward flow in the buffer layer whose magnitude increases with ζ. All these impact the near surface turbulence production processes as evidenced from the inner-scaled turbulence production, ⟨Pii⟩ + , which shows buffer layer modulation with increasing wave slope, ζ. In fact, we observe that ζ changes the proportion of form drag relative to viscous drag for a fixed roughness Reynolds number, a + .
The effect of surface undulations on the TBL is to enhance nearsurface mixing and reduce anisotropy of the buffer region at higher ζ. 35 Specifically, the surface undulations reduce the inner-scaled streamwise variance, ⟨u ′2 ⟩ + x,z,t , while enhancing the inner-scaled vertical (⟨v ′2 ⟩ + x,z,t ) and spanwise (⟨w ′2 ⟩ + x,z,t ) variances in the surface layer of the TBL. Thus, the flow becomes increasingly isotropic closer to the wall at higher steepness. For these shallow wavy surfaces, the streamwise variance is predominantly generated from shearinduced mechanisms, i.e., P u ′ v ′ 11 + driven by the strain rate term, d⟨u⟩ + z,t /dy + . Nevertheless, the surface-induced contribution, P u ′ u ′ 11 + , driven by d⟨u⟩ + z,t /dx + impacts the key aspects of the two-dimensional production structure, especially near the windward wavy surface. This secondary term can be related to production from surfaceinduced dispersive stresses and determines the key trends observed in the double-averaged (one-dimensional) statistics, ⟨P u ′ v ′ 11 ⟩ + x and ⟨u ′2 ⟩ + x,z,t . Spatially, the streamwise variance (⟨u ′2 ⟩ + z,t ) generation occurs in the leeward side of the wave in the buffer layer in response to the large inner-scaled gradients, d⟨u⟩ + z,t /dy + , in the wake of the wave crest with strong turbulent mixing. Along the windward region of the wavy surface, ⟨u ′2 ⟩ + z,t is produced close to the wall due to surface-slope-induced positive covariance, i.e., ⟨u ′ v ′ ⟩ + z,t > 0.
As the wave slope (ζ) increases, it enhances the strain rate terms d⟨u⟩ + z,t /dy + and d⟨u⟩ + z,t /dx + and thereby variance production. Unlike the TBL over a flat surface, small amounts of vertical variance are produced from Reynolds stress-mean strain rate interactions that arise purely from surface heterogeneity effects. In this case, the mean strain rate terms responsible for production, d⟨v⟩ + z,t /dy + and d⟨v⟩ + z,t /dx + , are both small with d⟨v⟩ + z,t /dy + ≫ d⟨v⟩ + z,t /dx + (since ζ ≪ 1). Therefore, as observed for ⟨u ′2 ⟩ + z,t , the vertical mean velocity gradient, d⟨v⟩ + z,t /dy + , determines the qualitative structure of the production term, P 22 + . Spatially, ⟨v ′2 ⟩ + z,t is produced along both the leeward (near the surface) and windward (near and away from the surface) sides of the wave with the latter being more dominant. Despite this production, the primary source of ⟨v ′2 ⟩ + z,t is through the pressure-rate-of-strain mechanism, which generates vertical fluctuations in the upper buffer layer. Near the surface, the same mechanism converts vertical fluctuations into streamwise and spanwise fluctuations on account of wall blockage. As the wave slope increases, both the double-averaged production, ⟨P 22 ⟩ + x , and pressure-rate-of-strain term, ⟨R 22 ⟩ + x , increase in magnitudes closer to the surface resulting in faster growth of ⟨v ′2 ⟩ + x,z,t (through the lower buffer layer) resulting in a downward shift in peak ⟨v ′2 ⟩ + x,z,t . As spanwise fluctuations are not as severely blocked by the wall as vertical fluctuations, ⟨w ′2 ⟩ + x,z,t grows faster (relative to vertical variance) with y + and thereby causes it to peak in the lower buffer layer. Due to spanwise homogeneity, the production of ⟨w ′2 ⟩ + z,t , P 33 + from Reynolds stress-strain rate interactions is zero. Therefore, generation of ⟨w ′2 ⟩ + z,t occurs through the pressure-rate-ofstrain term, R 33 + , which converts streamwise and vertical fluctuations into spanwise turbulence, especially along the windward side of the wavy surface. With an increase in the wave slope, this conversion process is enhanced, especially within the viscous and buffer layers. The 2D structure of the pressure-rate-of-strain terms near the surface shows three very distinct phase (ϕ)-dependent conversion mechanisms along the wave. These include splat-type phenomena along the leeward side and surface-induced generation of spanwise fluctuations over the windward side. Away from the surface, the well-known return-to-isotropy mechanism is observed in regions of strong vertical and streamwise variance. We also explore different metrics for quantifying the vertical extent of the surface (or roughness) layer, i.e., height beyond which surface effects are not observed. We look at the surface induced variance double-averaged production (⟨P u ′ u ′ 11 ⟩ + x , ⟨P 22 ⟩ + x ) as well as dispersion in this production structure (such as ⟨P u ′ u ′ 11 ⟩ + x − ⟨P u ′ u ′ 11 ⟩ + * ) to make these characterizations. The roughness layer height varies between ∼3a + and 5a + , which is larger than the earlier reported 61 values of ∼2a + using other metrics. More research is needed to assess such variability in estimates and Reynolds number sensitivity of these conclusions.