Reynolds number dependence of Reynolds and dispersive stresses in turbulent channel flow past irregular near-Gaussian roughness

Direct numerical simulations of fully-developed turbulent channel flow with irregular rough walls have been performed at four friction Reynolds numbers, namely, 180, 240, 360 and 540, yielding data in both the transitionally- and fully-rough regime. The same roughness topography, which was synthesised with an irregular, isotropic and near-Gaussian height distribution, is used in each simulation. Particular attention is directed towards the wall-normal variation of flow statistics in the near-roughness region and the fluid-occupied region beneath the crests, i.e. within the roughness canopy itself. The goal of this study is twofold. (i) Provide a detailed account of first- and second-order double-averaged velocity statistics (including profiles of mean velocity, dis- persive stresses, Reynolds stresses, shear stress gradients and an analysis of the mean force balance) with the overall aim of understanding the relative importance of “form-induced” and “turbulence-induced” quantities as a function of the friction Reynolds number. (ii) Investigate the possibility of predicting the levels of streamwise dispersive stress using a phenomenological closure model. Such an approach has been applied successfully in the context of idealised vegetation canopies (Moltchanov & Shavit, 2013, Water Resour. Res. , vol. 49, pp. 8222-8233) and is extended here, for the first time, to an irregular rough surface. Overall, the results reveal that strong levels of dispersive stress occur beneath the roughness crests and, for the highest friction Reynolds number considered in this study, show that the magnitude (and gradient) of these “form-induced” stresses exceed their Reynolds stress counterparts. In addition, this study emphasises that the dominant source of spatial heterogeneity within the irregular roughness canopy are “wake-occupied” regions and that a suitable parameterisation of the wake- occupied area is required to obtain an accurate prediction of streamwise dispersive stress.


Introduction
Practical fluid mechanics problems regularly involve turbulent flow past irregular rough surfaces. For instance, broadband roughness distributions affect the onset of laminar-turbulence transition on compressor blades (Goodhand et al., 2016), the pressure drop along steel pipes (Shockling et al., 2006) and the hydrodynamic resistance of marine vessels . In addition, irregular forms of surface roughness influence a wide range of physical processes in naturally-occurring turbulent flows. Some examples include: Mass transfer rates in coral reefs (Monismith, 2007), sediment transport in vegetated channels (Nepf, 2012) and turbulent transport mechanisms in gravel-bed rivers (Nikora and Smart, 1997). Considering their widespread ecological and technological importance, the fluid dynamic properties of irregular rough surfaces continue to be a rich area of research.
Over the past two decades, substantial effort has been dedicated to the study of turbulent flow past irregular rough walls. On the experimental side, hot-wire anemometry (HWA), particle image velocimetry (PIV) and laser doppler anemometry (LDA) have been employed to acquire turbulence measurements in irregular rough-wall pipes (Shockling et al., 2006;Langelandsvik et al., 2008), boundary layers Barros and Christensen, 2014;Morrill-Winter et al., 2017), channels (Flack et al., 2016;Barros et al., 2018) and openchannels (Nikora et al., 2001;Spiller et al., 2015;Nikora et al., 2019). Similarly, on the computational side, large-eddy simulations (LES) and direct numerical simulations (DNS) have been performed to study turbulent boundary layers (Anderson and Meneveau, 2011;Cardillo et al., 2013), channels (Napoli et al., 2008;De Marchis and Napoli, 2012;Busse et al., 2015;Thakkar et al., 2016) and open channels (Scotti, 2006;Yuan and Piomelli, 2014a;Forooghi et al., 2017) with irregular rough walls. Most of this past work has focussed on addressing three main issues. (i) Assessing the validity of Townsend's outer-layer similarity hypothesis (Townsend, 1976) with the aim of quantifying the interaction (or lack thereof) between the near-wall and outer flow. (ii) Identification of key topographical parameters, e.g. effective slope (Napoli et al., 2008) and skewness (Flack et al., 2016), with the aim of developing empirical methods to predict roughness effects in practical flows. (iii) Characterisation of turbulent secondary motions, e.g. the low-and high-momentum pathways observed by Barros and Christensen (2014), to quantify the degree of spatial heterogeneity on the cross-stream plane.
Several recent studies of turbulent flow past irregular roughness have directed attention towards the so-called "double-averaging" methodology and "form-induced" dispersive stresses. Double-averaging refers to a two-step method whereby consecutive temporal and spatial averages are applied to the Navier-Stokes equations. The first step yields the Reynolds-averaged Navier-Stokes (RANS) equations, whereas the second step yields the double-averaged (DA) Navier-Stokes equations. The Double-Averaged Navier-Stokes (DANS) equations were formulated in the context of atmospheric boundary layer studies to characterise turbulence above and within vegetation canopies (Raupach and Shaw, 1982). Detailed discussions regarding the DANS equations and dispersive stresses can be found in past work by Giménez-Curto and Lera (1996), Nikora et al. (2007) and Moltchanov et al. (2011). Dispersive stresses arise as a direct result of double averaging and represent momentum flux in the near-roughness region due to spatial heterogeneity in the time-averaged flow . The degree of heterogeneity is determined by the roughness height distribution, whose peaks induce low-momentum regions in their wake and zero velocity on their solid-fluid interface (Moltchanov and Shavit, 2013).
DA statistics (including both dispersive and spatially-averaged Reynolds stresses) have already received considerable attention in urban roughness studies (Cheng and Castro, 2002;Lien and Yee, 2004;Castro et al., 2006;Coceal et al., 2006;2007) and canopy flow studies (Poggi et al., 2004;Poggi and Katul, 2008;Moltchanov and Shavit, 2013), which tend to focus on regular rough surfaces comprised of geometric elements, e.g. cubes or cylinders. In the context of turbulent pipe flow, Chan et al. (2018) used DNS to compare the relative magnitude of dispersive and Reynolds stresses, noting that the former stresses scale on the wavelength of geometrically-scaled sinusoidal roughness topographies. The comparative lack of DA statistics in the context of irregular forms of surface roughness is conveyed in Table 1, where forty previous studies relevant to this work are listed. Brief details of the experimental methodology, flow configuration, roughness topography and the statistical quantities reported in each study are included for comparison. Some notable observations based on Table 1 include: (i) Single-point experimental techniques (e.g. HWA) are rarely used to accumulate DA statistics, possibly because of the time-intensive nature of obtaining spatially-resolved data. (ii) Whole flow-field experimental techniques (e.g. PIV) are more often used to calculate dispersive stresses, probably because the multi-dimensional nature of the instantaneous vector fields is better-suited to time-then-space averaging. (iii) With the exception of a number of recent LES and DNS campaigns, detailed comparisons of Reynolds and dispersive stresses (and their gradients) remain scarcely available. Dispersive stresses have been examined in the context of rough-wall turbulent open channel flow by Forooghi et al. (2018), whose DNS results show that the peak value of streamwise dispersive stress can become comparable to its Reynolds stress counterpart. A recent numerical experiment by Jelly and Busse (2018b) clarified the relative contributions of dispersive and Reynolds shear stress towards the rise in mean momentum deficit in a fully-developed turbulent channel flow with highly skewed irregular rough walls. In addition, Yuan and Jouybari (2018) performed DNS to compare Reynolds and dispersive stresses in a rough-wall turbulent open channel flow, highlighting how the magnitude of "forminduced" stress is strongly influenced by the spectral content of the roughness distribution. In each of these past studies, the Reynolds number was fixed and the effect of systematically varying the roughness topography upon the DA statistics was investigated.
The present study complements this past work by adopting the opposite approach: The roughness topography is fixed and relative magnitudes of the dispersive and Reynolds stress are investigated as a function of the Reynolds number. DNS of fully-developed turbulent channel flow with irregular rough walls have been performed at four separate friction Reynolds numbers. Particular attention is paid to the vertical variation of DA statistics in the near-roughness region and within the roughness canopy itself -regions where experimental data can be very challenging to obtain. The relative magnitudes of dispersive and Reynolds stresses (and their wall-normal gradients) are examined across a range of friction Reynolds numbers -enabling key trends to be identified and comparisons against past relevant works to be made. In addition, time-averaged stress distributions are also examined with the aim of complementing the DA data. Finally, this study takes the first steps towards modelling dispersive stresses in the context of irregular roughness by implementing and assessing the performance of a phenomenological closure model proposed in past work by Moltchanov and Shavit (2013).
This document is organised into four sections. Section 2 describes the computational aspects of this work including details of the surface generation algorithm, simulation setup and statistical averaging procedures. The results of this study are presented in Section 3 where the Reynolds number dependence of first-and second-order DA velocity statistics are examined and discussed. Finally, in Section 4, the conclusions of this work are given and recommendations for future research are made.

Computational aspects
This section describes the computational aspects of this work and is divided into three parts. First, the algorithm employed to synthesise an irregular surface with a near-Gaussian height distribution is described. Second, details of the computational setup and the key simulation parameters are provided. Lastly, the double-averaging methodology is introduced and the dispersive and Reynolds stress tensors are defined.

Surface characterisation
An irregular three-dimensional rough surface can be described using a two-dimensional heightmap, denoted here as h(x 1 , x 2 ), where x 1 and x 2 are the streamwise and spanwise directions, respectively. A heightmap quantifies the local elevation measured relative to a reference mean plane. Throughout this work the reference mean plane is located at = x 0, 3 where x 3 denotes the wall-normal direction. It is customary to characterise a height distribution using a set of amplitude parameters. Four commonly used amplitude parameters include the mean absolute height, Sa, and the root-mean-squared height, Sq, defined here as as well as skewness, Ssk, and kurtosis, Sku, defined here as where A is the planform area of the heightmap. An extensive discussion of these four parameters (and many more) can be found in past work by Mainsah et al. (2001) and Leach (2013). For any naturally-occurring or synthetic rough surface, all possible combinations of skewness and kurtosis are bounded by the Pearson (1916) inequality Jelly and A. Busse International Journal of Heat and Fluid Flow 80 (2019) The quadratic boundary defined by inequality 3 is plotted in Fig. 1, along with one-hundred and seventy-four skewness-kurtosis combinations taken from eleven independent surface measurement campaigns. As expected, each measurement point falls within the permissible bounds of the Pearson (1916) inequality 3. Furthermore, the majority of the measured data is clustered around the Gaussian case, corresponding to zero skewness = Ssk ( 0) and kurtosis equal to three = Sku ( 3). The remaining data points include several notable outliers, some with extremely non-Gaussian properties, which, based on the available data, mostly lie in the lower half-plane (Ssk < 0). Nonetheless, the majority of the measured data points are near-Gaussian. Motivated by this observation, this work focuses on the fluid dynamic properties of a "custom-made" surface synthesised with a normal height distribution.

Surface generation algorithm
An irregular surface with was generated by taking weighted linear combinations of a Gaussian random number matrix using a moving average process. The surface was synthesised with an isotropic exponential autocorrelation coefficient function and doubly-periodic boundaries. The same surface generation algorithm was used in recent studies by Jelly and Busse (2018a,b) and is based on past work by Patir (1978). A smoothly varying heightmap was obtained by applying a two-dimensional low-pass Fourier filter to the discrete point cloud using the method described by Busse et al. (2015) -the resulting filtered surface is shown in Fig. 2. After filtering, the heightmap was scaled to have a mean-peak-to-valley height of where δ is the mean channel half height. To compute this quantity, the heightmap was split into 5 × 5 equally-sized tiles and the maximum and minimum height of each tile was evaluated. S z,5 × 5 is defined as the difference between the mean of the maxima and mean of the minimafurther details can be found in Thakkar et al. (2016).
The key topographical parameters of the filtered surface are provided in Table 2. As was previously mentioned, the surface has near-Gaussian properties, i.e. negligible skewness (Ssk ≈ 0) and kurtosis equal to three = Sku ( 3), which is representative of what is often encountered in reality (see Fig. 1). In addition to the amplitude parameters listed in Table 2, the streamwise effective slope, ES, defined as is also included. Past work by Napoli et al. (2008) and Schultz and Flack (2009) Table 2, the surface in this study has an effective slope of = ES 0.37, and, as a result, the roughness effects reported herein fall into the so-called "roughness flow regime".

Simulation set up
DNS of rough-walled incompressible turbulent channel flow were performed using an iterative embedded-boundary algorithm (Busse et al., 2015). The current DNS algorithm is thoroughly validated and has been used in past studies related to the this work, e.g. Thakkar et al. (2018) and Jelly and Busse (2018a,b). A schematic of the rough-wall channel flow configuration adopted throughout this study is shown in Fig. 3.
The velocity components in the streamwise (x 1 ), spanwise (x 2 ) and wall-normal (x 3 ) directions are u 1 , u 2 , and u 3 , respectively, and p is the fluctuating pressure. The irregular surface shown in Fig. 2 was used to enforce impermeable, no-slip boundary conditions on both the upper and lower walls. Periodic boundaries were applied in the streamwise and spanwise directions. A constant (negative) mean streamwise pressure gradient, Π, was prescribed to drive the flow through the channel. The mean streamwise pressure gradient and mean friction velocity, u τ , are related through the formula where δ is the mean channel half-height and ρ is density. The friction Reynolds number is defined as Re τ ≡ u τ δ/ν, where ν is kinematic viscosity. Throughout this document, inner-scaled quantities are marked with a superscript +. The inner-scaled wall-normal position is measured relative to the mean plane of the lower surface = x ( / 0) 3 and is denoted as . Four rough-wall simulations were conducted at four different friction Reynolds numbers, namely, 180, 240, 360 and 540. Additional smooth-wall calculations were performed at each friction Reynolds number -giving a total of eight separate simulations. Key computational parameters for the rough-and smooth-wall cases are provided in Table 3. The current mesh resolution for the rough-wall simulations adheres to the guidelines and recommendations made by Busse et al. (2015). For all rough-wall simulations considered in this study, statistical data was accumulated for a minimum of seventy-five      Tennekes and Lumley (1972), the near-wall layer over a rough surface can be described by two length-scales: (i) A characteristic roughness height, k, and (ii) the viscous length-scale, defined here as ℓ ν ≡ ν/u τ . If the ratio of (i) and (ii) is defined as the roughness Reynolds number, i.e. + k ku / , then roughness effects will commence as the viscous length scale and the characteristic roughness scale become commensurate, i.e. + k (1). Whereas the viscous length-scale is explicitly defined by the fluid properties and flow conditions, multiple definitions of k can exist for surfaces with a broad range of topographical scales. The inner-scaled roughness parameters corresponding to the irregular heightmap shown in Fig. 2 are provided in Table 4, along with the ratio of the mean channel-half height to the outer-scaled ,5 5 . The maximum and minimum values of k i differ by a factor of seven -highlighting the difficulty in selecting the "correct" k to characterise the irregular near-Gaussian surface under consideration in this work. For example, considering that the inner-scaled mean absolute height, + Sa , is commensurate with the thickness of a smooth-wall viscous sublayer for the Re τ range of this study, then a very weak roughness effect would be expected. On the other hand, the inner-scaled mean-peak-to-valley height, × + S , z,5 5 is a factor of seven greater than + Sa and implies a far stronger roughness effect should be anticipated. This ambiguity also applies to the ratio of the mean channel half-height to the outer-scaled roughness parameters, δ/k i . For roughness with δ/k < 40, Jiménez (2004) argues that the defining mechanisms of near-wall turbulence no longer exist, ultimately leading to the breakdown of wall similarity. However, this limit was made in the context of two-dimensional transverse square bar roughness, where k is defined explicitly by the bar height. In contrast, there is currently no universally-accepted definition of k for irregular three-dimensional rough surfaces. As a result, defining an unambiguous threshold on δ/k for a broadband roughness distribution is not straightforward. Nonetheless, the range of δ/k i corresponding to the irregular surface under consideration in this study (see Table 4) is comparable to past relevant work. For instance, Thakkar et al. (2016) reports the ratio of the mean absolute height to the mean-peak-to-valley height spans the range 4 ≲ Sa/S z,5 × 5 ≲ 7 for a set of seventeen surface scans -comparing well against the values listed in Table 4.
Bulk flow properties corresponding to the smooth-and rough-wall simulations are compared in Table 5. Since the mean streamwise pressure gradient was held constant throughout each simulation, an increase in drag can be inferred from a decrease in centreline Reynolds number, Re cl ≡ U cl δ/ν, bulk Reynolds number, Re b ≡ U b δ/ν, or, alternatively, an increase in the roughness function, + U . In this work, the roughness function is defined as the difference between the smoothand rough-wall centre-line streamwise velocity at a matched friction Reynolds number. The roughness function is plotted against its equivalent sandgrain roughness, + k , s in Fig. 4. The equivalent sandgrain roughness was obtained by collapsing the value of 540 on to the fully-rough asymptote. Fitting the data in this way reveals that the equivalent sandgrain roughness and the mean peak-to-valley height are related through the formula = + × + k S 1.06 s z,5 5 . This relationship shows that × + S z,5 5 is a better candidate to estimate + k s of an irregular Gaussian surface relative to other available roughness length-scales, e.g. the mean absolute height which is a factor of seven smaller (see Table 4). The approximation + × + k S s z,5 5 was also noted in past work by Busse et al. (2017) and Thakkar et al. (2018), who examined the Reynolds number dependence of mildly non-Gaussian irregular roughness. The onset of the fully-rough regime is widely accepted to occur at an equivalent sandgrain roughness of + k 70 s (Flack and Schultz, 2010), corresponding to a roughness function of + U 7 based on Nikuradse's sand-grain data (Nikuradse, 1933). Based on this criterion, the highest friction Reynolds number case considered in this study can be classified as fully rough, whereas all other cases fall into the transitionally rough regime. Whilst some basic performance trends can be drawn from Table 5 and Fig. 4, better understanding can be obtained by evaluating and examining first-and second-order double-averaged (DA) velocity statistics.

Double-averaging methodology
Throughout this work, statistical quantities were computed using a double-averaging approach (Raupach and Shaw, 1982). Any instantaneous field variable, say θ, can be DA by the successive application of the time-averaging operator and the spatial-averaging operator In Eq. (7), the void fraction function, represents the ratio of the fluid-occupied area, A f , to the total area of the wall-parallel plane, = A L L 1 2 . In the region extending below the highest roughness crest (x 3 /h max < 1) the void fraction function is less than one. In solidoccupied regions instantaneous field variables are set equal to zero, i.e.
= t x ( , ) 0. Averaging the flow in this way ensures that only fluid-occupied points contribute towards DA statistics beneath the roughness crests -the so-called "intrinsic average" (Gray and Lee, 1977). Alternatively, if no distinction between fluid-and solid-occupied points is made beneath the crests, then spatially averaged data in this region corresponds a "superficial average", which, with reference to Eq. (7), can be obtained by prescribing a unit void fraction function, i.e. = 1. The former approach is adopted throughout this study and, hence, all DA data presented herein corresponds to the intrinsic average.
As noted by Nikora et al. (2007), the void fraction function, ϕ, is equivalent to the cumulative density function (CDF) of the heightmap, which, for a standard Gaussian distribution, can be expressed as where erf denotes the error function. Wall-normal profiles of the standard Gaussian CDF (Eq. (8)) and the void fraction function are compared in Fig. 5. Overall, the two profiles match closely at all wallnormal positions and emphasise the smoothly-varying, near-Gaussian properties of the surface synthesised to undertake this study. Note that  Table 5. The Colebrook formula, ) and Nikuradse's fully-rough asymptote, ) are also shown for reference, both of which were computed using a von Kármán coefficient of = 0.41. Nikuradse (1933) sand-grain data ( ) is also included.

T.O. Jelly and A. Busse
International Journal of Heat and Fluid Flow 80 (2019) 108485 for some forms of regular roughness, e.g. transverse square bars (Krogstad et al., 2005) or vertical cylinder arrays (Moltchanov and Shavit, 2013), the void fraction function becomes constant beneath the crests and exhibits a discontinuity at the canopy edge due to a sudden change in the fluid-occupied area. As noted by Pokrajac and Manes (2008), a discontinuous void fraction function can complicate the spatial differentiation of DA data. However, as shown in Fig. 5, the irregular surface in this work has a smooth void fraction function, and, as a result, spatial gradients of DA statistics can be calculated using standard differencing techniques.
Considering the time-averaging operator 6 and the spatial-averaging operator 7, the DA methodology leads to the following triple decomposition of an instantaneous field variable where the "form-induced" dispersive component, x ( ), is defined as the difference between the local time-average and the DA, i.e.
x x x ( ) ( ) ( ), 3 and θ′(x, t) denotes the turbulent fluctuation. Based on the triple decomposition 9, the local Reynolds stress tensor is defined here as The spatial-average of the Reynolds stress tensor is defined here as Finally, the dispersive stress tensor is defined here as

Results
This section presents the key results of this study and is divided into five parts. First, inner-and outer-scaled profiles of DA axial velocity are presented and examined, along with an analysis of the mean reverse flow that occurs within the roughness canopy. Second, the relative magnitude of dispersive and Reynolds normal stresses are compared in the near-roughness region and in the region extending below the crests. Third, the relative magnitude of dispersive and Reynolds shear stresses (and their gradient) are evaluated and compared, along with a breakdown of the mean force balance above the roughness crests. Fourth, time-averaged contours of dispersive and Reynolds stress within the roughness canopy are examined in order to complement the analysis of DA data. Lastly, a closure model for the streamwise component of dispersive stress -devised first by Moltchanov and Shavit (2013) -is implemented and its accuracy and robustness are quantified.

Analysis of DA axial velocity profiles
Inner-scaled profiles of DA streamwise velocity normalised by the mean friction velocity are compared at each friction Reynolds number in Fig. 6. For the rough-wall data, the log-law emerges in the immediate vicinity of the viscous-scaled crests + + x h ( ) 3 max and then follows the rough-wall log-law where κ is the von Kármán constant and B is the smooth-wall intercept. The preservation of the log-law is in-line with Townsend's outer-layer similarity hypothesis (Townsend, 1976), and is consistent with the notion that the mean flow in the outer region becomes independent of the viscous-scaled surface condition. With reference to Table 4, the height of the highest roughness crest is approximately equal to four times the standard deviation of the heightmap, i.e. h max /Sq ≈ 4. This factor of four agrees well with past work that states roughness effects typically extend up to a "few roughness heights" from the wall (Raupach et al., 1991;Nikora et al., 2001;Flack et al., 2007). Relative to the smooth-wall data, the downward shift of the rough-wall log-law increases with increasing friction Reynolds number -the corresponding values of the roughness function, + U , are given in Table 5.   Table 3. All data has been normalised using the friction velocity, u τ . The height of viscous-scaled roughness crest for each rough-wall case is denoted using the vertical gray line. The inset panel shows the wall-normal variation of the roughness function.

Table 3
Simulation parameters for smooth-and rough-wall cases including: Friction Reynolds number (Re τ ); domain length (L 1 ); domain width (L 2 ); viscous-scaled mesh spacings including streamwise The wall-normal variation of the roughness function is plotted in the inset panel of Fig. 6. In the outer region, the roughness function achieves a constant value -a necessary condition for outer-layer similarity and, hence, the preservation of the log-law (Eq. (13)). It is worth noting that the roughness function achieves its maximum value below the height of the highest roughness crest (x 3 /h max < 1) before decaying towards zero as the mean roughness plane = x ( / 0) 3 is approached. A comparable overshoot of the roughness function is evident in past work by Chan et al. (2015) who simulated turbulent flow through pipes with a three-dimensional sinusoidal roughness topography. In contrast, such an overshoot is not guaranteed for surfaces with high levels of negative skewness (e.g. see "pits-only" surface studied by Jelly and Busse, 2018b) where + U decays monotonically with decreasing wall-normal position. This observation implies that the roughness peaks are responsible for the increased mean momentum deficit below the crests.
Outer-scaled profiles of DA streamwise velocity normalised by the mean friction velocity are compared at each friction Reynolds number in Fig. 7. Relative to the smooth-wall data, the rough-wall data show a clear reduction in centre-line and bulk velocities at each friction Reynolds number -consistent with the bulk flow quantities listed in Table 5. A close-up view of the reverse flow region that develops within the roughness canopy is shown on the inset panel in Fig. 7. The peak magnitude of reverse flow increases as friction Reynolds number increases and its position moves from right to left, i.e. downwards, deeper within the cavities. For example, when the friction Reynolds number is increased from 180 to 540, the absolute peak value of reverse flow increases by a factor of five and approaches 2.5% of the centre-line velocity. These observations agree with the past findings of Busse et al. (2017), who noted that high-amplitude reverse flow events are more likely to occur in the bottom half of irregular roughness canopies for increasing friction Reynolds numbers on the range 90 < Re τ < 720.
To further characterise the reverse flow region within the roughness canopy, isosurfaces of negative time-averaged streamwise velocity are compared at = Re 180 and 540 in Fig. 8. The isosurfaces are coloured by wall-normal position and are shown on the same (2 × 2)/δ sub-domain. Comparing the higher friction Reynolds number data against that of the lower, the localised pockets of reverse flow appear to spread laterally (along the x 2 direction) and coalesce to form elongated strips. A larger volume of mean recirculation therefore exists at higher values of Re τ , and, as a result, a stronger DA reverse flow is established (see inset panel on Fig. 7).
The preceding analysis of first-order axial velocity statistics can be summarised as follows. (i) Outer similarity is achieved in the immediate vicinity of the highest roughness crest. (ii) As a consequence of (i), the log-law emerges just above the viscous-scaled crests and shifts downwards with increasing friction Reynolds number (Fig. 6). (iii) Below the roughness crests, a DA reverse flow region occurs within the bottom half of the roughness canopy whose size and intensity both increase as the friction Reynolds number becomes larger (Figs. 7 and 8). Next, the Reynolds-number dependence of second-order velocity statistics will be examined in the context of the dispersive and spatially-averaged Reynolds normal stresses.

Analysis of dispersive and spatially-averaged Reynolds normal stresses
As was previously mentioned, the relative magnitude of dispersive and Reynolds stresses has been reported in recent work relevant to this study, e.g. Forooghi et al. (2018) and Yuan and Jouybari (2018). In those works, the Reynolds number was held fixed and the effect of varying the roughness topography was investigated. In this work, the relative magnitude of dispersive and Reynolds stresses are compared at four friction Reynolds numbers for the same roughness topography.
Line profiles of the dispersive and spatially-averaged Reynolds normal stresses for each friction Reynolds number are shown in Fig. 9.
The left-hand column shows rough-wall data on the range in order to emphasise the relative magnitude of dispersive and Reynolds stresses in the near-roughness region and within the roughness canopy itself. The inset panels show the peak stress magnitude as a function of the friction Reynolds number, where smooth-wall data points are included for reference. The right-hand column shows profiles of the smooth-and rough-wall Reynolds stresses plotted against inner-scaled wall normal position, + x , 3 in log-linear format. All normal stresses have been normalised using the square of the mean friction velocity, u 2 .
Profiles of streamwise dispersive stress, D 11 , and spatially-averaged Reynolds stress, < R 11 > , are compared in Fig. 9(a). The first notable observation based on this data is the collapse of the Reynolds stress profiles above a height of x 3 /h max ≳ 4. The profiles of smooth-and rough-wall Reynolds stress also collapse onto a single curve in the same   Fig. 7. Outer-scaled DA streamwise velocity profiles. Line-types are specified in Table 3. All data has been normalised using the friction velocity, u τ . The height of the highest roughness crest is denoted as the vertical gray line. The inset panel shows the mean recirculation region within the roughness canopy.
T.O. Jelly and A. Busse International Journal of Heat and Fluid Flow 80 (2019) 108485 region for each value of Re τ ( Fig. 9(b)). These observations are in-line with Townsend's outer-layer similarity hypothesis (Townsend, 1976) and complement the collapse of first-order statistics discussed in the previous subsection (see Fig. 6). A second notable observation based on Fig. 9(a) and (b) is how the peak value of < R 11 > varies with the friction Reynolds number for the smooth-and rough-wall data. Whereas the peak value of < R 11 > for the smooth-wall data shows a weak dependence on the friction Reynolds number (for a detailed discussion see work by Lee and Moser, 2015), the rough-wall data shows a far stronger trend. For example, increasing the friction Reynolds number from 180 to 540 results reduces the peak value of < R 11 > by one third, which is related to the disruption of near-wall streaks and quasistreamwise vortices in the vicinity of the roughness crests (Krogstad et al., 2005;Schultz and Flack, 2007). In addition, the position of peak < R 11 > moves from right to left and occurs beneath the roughness crests above a friction Reynolds number of 360. Relative to the Reynolds stress, the levels of dispersive stress remain negligible at all heights above the roughness crests (x 3 /h max > 1). However, in the region just below the crests, the dispersive stress begins to increase and the profiles of D 11 collapse onto a single curve on the range 0.6 ≲ x 3 / h max ≲ 1. A similar collapse of streamwise dispersive stress profiles has been reported by Chan et al. (2018), who simulated turbulent pipe flow at a friction Reynolds number of 540 with geometrically-scaled sinusoidal rough walls. Furthermore, at a friction Reynolds number of 540, the peak value of D 11 exceeds its Reynolds stress counterpart on the range 0 ≲ x 3 /h max ≲ 0.5 -underlining the strong deviations from streamwise-spanwise homogeneity within the roughness canopy. Finally, with reference to the inset panel of Fig. 9(a), we note that the peak magnitudes of dispersive and Reynolds stress converge towards a common value as the friction Reynolds number increases. As will be shown later, the primary source of streamwise dispersive stress are lowmomentum regions that develop in the immediate downstream vicinity of the protruding roughness crests, i.e. "wake-occupied" regions. Profiles of spanwise dispersive stress, D 22 , and spatially-averaged Reynolds stress, < R 22 > , are compared in Fig. 9(c). In contrast to the streamwise Reynolds stresses ( Fig. 9(a)), the levels of < R 22 > intensify for increasingly large friction Reynolds numbers in the vicinity of the roughness crests. The smooth-wall profiles show the same trend ( Fig. 9(d)), and, when plotted against inner-scaled wall-normal position, collapse on top of the rough-wall data in the outer-region. Furthermore, the peak value of < R 22 > for the smooth-and rough-wall data agree to within 5% for each value of Re τ (see inset panel on Fig. 9(c)). A comparable match has been reported by Yuan and Piomelli (2014a) in LES of smooth-and a rough-wall open-channel flows at a friction Reynolds numbers of 180 and 1000. Similar to the streamwise dispersive stress (Fig. 9a), non-negligible levels of D 22 only exist beneath the roughness crests and begin to exceed the local value of < R 22 > at friction Reynolds number of 540. Finally, whereas the peak value of D 22 increases monotonically with increasing friction Reynolds number, it remains a factor of two smaller than its "turbulence-induced" counterpart for each value of Re τ .
Profiles of wall-normal dispersive stress, D 33 , and spatially-averaged Reynolds stress, < R 33 > , are compared in Fig. 9(e). Overall, the data share several similarities with the spanwise stresses shown in Fig. 9(c) and (d). These include negligible levels of D 33 above the roughness crests, non-negligible levels of D 33 within the roughness canopy, relative to < R 33 > , and a monotonic increase in the peak value of D 33 with increasing friction Reynolds number. Relative to the streamwise and spanwise dispersive stresses, however, the magnitude of D 33 is notably smaller -which in is agreement with the past results of Busse et al. (2015) and Forooghi et al. (2018). Profiles of smooth-and rough-wall < R 33 > plotted in Fig. 9(f) show an excellent level of collapse in the outer region, and become indistinguishable above the viscous-scaled crests > A similarly striking collapse in wallnormal Reynolds stress profiles has been reported by Chan et al. (2018) in turbulent pipe flow with three-dimensional sinusoidal rough walls. On the other hand, Volino et al. (2011) noted that < R 33 > became 20% higher in the outer region of a turbulent boundary-layer past "large" two-dimensional bar roughness with δ/k ≈ 33, compared to the smooth-wall case. However, as demonstrated in later work by Krogstad and Efros (2012), outer layer similarity on < R 33 > was recovered for two-dimensional rough walls at higher Reynolds numbers and higher values of δ/k.
Three main conclusions can be drawn from the stress profiles plotted in Fig. 9: (i) Increasing the friction Reynolds number leads to an appreciable suppression of < R 11 > and a simultaneous enhancement of < R 22 > and < R 33 > in the vicinity of the roughness crests. This behaviour is consistent with the notion that turbulent flow in the nearroughness region is driven towards a more isotropic state, relative to the smooth-wall state (see Smalley et al., 2002 for further details). (ii) Non-negligible levels of dispersive stress occur beneath the roughness crests, i.e. within the roughness canopy itself. This behaviour agrees with a past urban roughness study by Cheng and Castro (2002) who noted that, compared to spatially-averaged Reynolds stresses, dispersive stresses become negligible above cubical obstacles. (iii) When plotted against inner-scaled wall-normal position, the Reynolds normal stresses exhibit excellent levels of outer-layer similarity at each friction Reynolds number. This observation is in good agreement with findings in recent numerical studies relevant to this work (Chan et al., 2018;Forooghi et al., 2018;Yuan and Jouybari, 2018), as well as a number of experimental campaigns (Wu and Christensen, 2007;Squire et al., 2016;Pathikonda and Christensen, 2017).

Analysis of dispersive and spatially-averaged Reynolds shear stresses
To complement the analysis of the normal stresses, dispersive and spatially-averaged Reynolds shear stress profiles are plotted in Fig. 10. The data are presented in the same format as Fig. 9, where the left-hand panel shows dispersive and Reynolds shear stresses on the range < < x h 1 / 5 3 max and the right-hand panel shows smooth-and rough- . All data has been normalised using the friction velocity, u τ . Note that only a (2 × 2)/δ sub-section of the full surface is shown.

T.O. Jelly and A. Busse
International Journal of Heat and Fluid Flow 80 (2019) 108485 Fig. 9. Profiles of dispersive and spatially-averaged Reynolds normal stresses. Left-hand column shows rough-wall data including spatially-averaged Reynolds normal stresses (black lines), and dispersive normal stresses (blue lines) including (a) streamwise, < R 11 > and D 11 ; (c) spanwise, < R 22 > and D 22 and (e) wall-normal, < R 33 > and D 33 , components. Inset panels on (a), (c) and (e) show maximum value of dispersive stress ( ), rough-wall Reynolds stress (∘) and reference smooth-wall data ( ). Right-hand column shows same rough-wall data as the left-hand column (black lines), along with smooth-wall Reynolds normal stresses (gray lines) including (b) streamwise, < R 11 > ; (d) spanwise, < R 22 > and (f) wall-normal, < R 33 > , components. Line-types correspond to different values of Re τ and are specified in Table 3. All data has been normalised using the friction velocity, u τ . The heights of the highest roughness crests are shown as the vertical gray lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

T.O. Jelly and A. Busse
International Journal of Heat and Fluid Flow 80 (2019) 108485 wall Reynolds stresses plotted against + x 3 in log-linear format. Again, the shear stress profiles have been normalised using the square of the mean friction velocity, u 2 .
Profiles of the dispersive and Reynolds shear stress are compared in Fig. 10(a). Beyond a height of x 3 /h max > 4, the Reynolds shear stress profiles collapse onto a straight line -indicating that the flow is converged for each value of Re τ considered in this study. Similar to the dispersive normal stresses ( Fig. 9(a),(c) and (e)), significant levels of dispersive shear stress only occur within the roughness canopy. In comparison to < R 13 > , the peak value of D 13 is at least a factor of four smaller for all friction Reynolds numbers considered here -an observation that suggests "form-induced" shear stress has a relatively weak influence upon the flow. With reference to the inset panel on 10 a, the peak value of rough-wall Reynolds shear stress matches closely with that of the smooth-wall and agrees to within 7% for all values of Re τ . In addition, the rough-wall Reynolds shear stress profiles collapse onto the smooth-wall data in the immediate vicinity above the viscous-scaled crests > + + x h ( ) 3 max and match very closely all the way to the channel centre ( Fig. 10(b)).
As noted by Manes et al. (2008), the magnitude of spatially-averaged Reynolds shear stress (or dispersive shear stress) is not enough to come to conclusion about their physical significance, since it is their wall-normal gradients that appear in the mean force balance. For the current channel-flow configuration, the DA streamwise momentum equation above the highest roughness crest (x 3 /h max > 1) can be written as where the bracketed term on the right-hand represents the combined effect of dispersive and spatially-averaged Reynolds shear stress gradients. Below the roughness crests (x 3 /h max < 1), extra terms enter the DANS equations related to additional viscous and pressure forcesfurther details can be found in Raupach and Shaw (1982), Nikora et al. (2007) or Moltchanov et al. (2015). The DA force balance above the roughness crests for a friction Reynolds number of 540 is shown in Fig. 11, where each term on the right-hand side of Eq. (14) is plotted for comparison. No data is shown beneath the roughness crests since Eq. (14) is invalid in this region. Above a height of approximately three roughness crests (x 3 /h max ≳ 3), the mean streamwise pressure gradient is balanced by the action of the Reynolds shear stress gradient, since both the viscous force and dispersive shear stress gradient become negligible in the outer region. As a result, the mean force balance Eq. (14) reduces to = in the outer flow (as for the smooth-wall case). As the roughness crests is approached, the Reynolds shear stress gradient goes through a zero- . Inset panel shows maximum value of dispersive shear stress ( ), rough-wall Reynolds shear stress (∘) and reference smooth-wall data ( ). (b) Inner-scaled profiles of rough-wall Reynolds shear stress (black lines) and smooth-wall Reynolds shear stress (gray lines). Line-types correspond to different values of Re τ and are specified in Table 3. All data has been normalised using the friction velocity, u τ . The heights of the highest roughness crests are shown as the vertical gray lines.. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) ) and Reynolds shear stress gradient, The sum of the terms (with the correct signs) is shown as open circles (∘). All data has been normalised using the friction velocity, u τ , and the height of the highest roughness crest, h max .
T.O. Jelly and A. Busse International Journal of Heat and Fluid Flow 80 (2019) 108485 crossing (switching its sign from negative to positive) and serves as a momentum source beneath this point. Furthermore, at the zero-crossing point of the Reynolds shear stress gradient, the dispersive shear stress gradient is the only remaining term, along with the viscous force, that maintains the flow against the mean pressure gradient -underlining the critical role that "form-induced" momentum transport plays in the near-roughness region. In the immediate vicinity of the roughness crests (x 3 /h max ≈ 1), the dispersive and Reynolds shear stress gradients are both order one quantities -implying their ratio is also of unit order. In contrast, the ratio of their magnitudes in this region is ten-to-one in favour of the Reynolds shear stress (see Fig. 10(a)) -emphasising the importance of considering shear stress gradients, as opposed to their magnitudes, when considering their relative contribution to the mean dynamics of the near-wall flow. Profiles of the dispersive and spatially-averaged Reynolds shear stress differentiated with respect to wall-normal height, x 3 , are compared at each friction Reynolds number in Fig. 12. As shown previously in Fig. 11, the gradient of dispersive shear stress becomes negligible above a height of approximately three roughness crests (x 3 /h max ≳ 3). Whereas the Reynolds shear stress gradient remains positive beneath the roughness crests, the dispersive shear stress gradient switches sign (from negative to positive), a behaviour that has also been observed in past DNS simulations of turbulent flow past urban-like roughness (Coceal et al., 2006) and open-channel experiments with periodically arranged cubes (Florens et al., 2013). The inset panel in Fig. 12 shows how the peak magnitude of dispersive and Reynolds shear stress gradient varies with the friction Reynolds number, where smooth-wall data is included for comparison. As the friction Reynolds number increases, the relative difference in the peak magnitude of the dispersive and Reynolds shear stress gradient decreases. This difference is notably smaller than the relative difference of the shear stress magnitudes (see inset panel in Fig. 10). For example, at = Re 540, the peak magnitude of dispersive and Reynolds shear stress differs by 73%, whereas the peak magnitude of their gradient differs by 34%. Extrapolation of this data suggests that the peak magnitude of the dispersive and spatiallyaveraged Reynolds shear stress gradient will be approximately equal at a friction Reynolds number of Re τ ≈ 1200.

Spatial heterogeneity of time-averaged quantities within the roughness canopy
To complement the preceding analysis of spatially-averaged velocity statistics (Figs. 9-12), the spatial distribution of time-averaged quantities in the fluid-occupied region below the roughness crests can also be examined. Analysis of such data can help to clarify how individual roughness elements influence of the local level of spatial heterogeneity in the time-averaged field, which, in turn, can provide a better understanding of the Reynolds-number trends drawn from spatially-averaged data.
Wall-normal slices of the streamwise Reynolds stress, = R u u , 11 1 1 and streamwise dispersive flux, u ũ˜, 1 1 are shown in Fig. 13(a),(b),(e) and (f). Note that the spatial average of u ũ1 1 is equal to the streamwise dispersive stress, i.e. = D u ũ1 1 1 1 (Eq. (12)). The data on the top row corresponds to a friction Reynolds number of 180, whereas the bottom row corresponds to a friction Reynolds of 540. The data on both rows correspond to a wall-normal position of = x h / 0.47 3 max -a height where the line profiles of D 11 and < R 11 > agree to within 20% (see Fig. 9(a)). Despite having comparable spatially-averaged values, the spatial distribution of u ũ1 1 and R 11 show a number of remarkable differences. First, whereas the highest levels of streamwise dispersive flux are concentrated within the wakes of the protruding roughness elements, the same regions correspond to Reynolds stress "dead-zones". Likewise, whereas the highest levels of Reynolds stress occur in the spaces between the roughness peaks, the same regions correspond to dispersive flux "dead-zones". Put in other words, regions of strong streamwise dispersive flux correlate with regions of weak streamwise Reynolds stress. As the friction Reynolds number is increased from 180 to 540, the intensity of u ũ1 1 in the wakes increases, which, ultimately, leads to stronger streamwise dispersive stress within the roughness canopy (see Fig. 9(a)). Increasing the friction Reynolds number also leads to a considerable reduction in the area occupied by Reynolds stress "dead-zones", which, again, explains the increased levels of < R 11 > in this region. In addition, the contours plotted in Fig. 13 show that the roughness peaks induce spatial heterogeneity within the canopy and highlight how the averaged effect of "wake-occupied" regions is the dominant source of < D 11 > .
Wall-normal slices of Reynolds shear stress, = R u u , 13 1 3 and streamwise-wall-normal dispersive flux, u ũ˜, 1 3 are compared in Fig. 13(c),(d),(h) and (g). The spatial distribution of u ũ1 3 reveals a scattered pattern of positive, negative and near-zero values that is, generally speaking, quite complicated. The irregular nature of u ũ1 3 is consistent with past work by Pokrajac et al. (2007), whose quadrant map analysis demonstrated that the spatial distribution of dispersive flux is determined by the nature of the underlying roughness topography. Unlike the contours of u ũ1 1 (Fig. 13(a) and (d)), there is no obvious connection between "wake-occupied" regions and the generation of u ũ1 3 . On the other hand, the contours of Reynolds shear stress show "dead zones" clustered around the roughness peaks ( Fig. 13(d) and (h)). In a manner similar to the streamwise Reynolds stress contours, the Reynolds shear stress "dead zones" shrink in size as the friction Reynolds number increases from 180 to 540. In addition, for the majority of the fluid-occupied region, the Reynolds shear stress remains negatively correlated in the spaces between the roughness elements. However, localised patches of positive Reynolds shear stress, i.e. R 13 > 0, are evident on the windward faces of the roughness peaks. Similar observations were made by De Marchis et al. (2010), who noted localised regions of positive Reynolds shear stress on the up-slope side of irregular two-dimensional surface roughness. However, the integral contribution of regions where R 13 > 0 towards the DA value is small, and, as a result, < R 13 > remains negative. Overall, the time-averaged stresses plotted in Fig. 13 show that the levels of spatial heterogeneity within the roughness canopy exhibit a strong Reynolds-number dependence, which, in turn, influences the levels of dispersive and spatially-averaged Reynolds stresses (Figs. 9 and 10). In addition, the Fig. 12. Profiles of dispersive and spatially-averaged Reynolds shear stress gradient. Outer-scaled profiles of the Reynolds stress gradient, < ∂R 13 /∂x 3 > (black lines), and dispersive stress gradients, < ∂D 13 /∂x 3 > (blue lines). Inset panel shows maximum value of dispersive shear stress gradient ( ), rough-wall Reynolds shear stress gradient (∘) and reference smooth-wall data ( ). The constant (negative) streamwise mean pressure gradient = ( 1) is shown as the horizontal gray line. All data has been normalised using the friction velocity, u τ . The height of the highest roughness crest is shown as a vertical gray line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

T.O. Jelly and A. Busse
International Journal of Heat and Fluid Flow 80 (2019) 108485 contour plots reveal that "wake-occupied" regions play a key role in determining the local levels of both u ũ1 1 and D 11 .

Assessment of Moltchanov and Shavit (2013) model
A previous experimental campaign by Moltchanov and Shavit (2013) demonstrated that "wake-occupied" regions are the dominant source of steamwise dispersive stress above (and within) idealised vegetation canopies. In that work, it was argued that if the time-averaged streamwise velocity becomes sufficiently small in the immediate vicinity downstream of the roughness elements, i.e. u 0, 1 then the magnitudes of dispersive and DA velocity in a "wake- (b,f) streamwise Reynolds stress, + R 11 ; (c,g) streamwise-wall-normal dispersive flux, + + u ũ1 3 and (d,h) Reynolds shear stress, + R 13 ; The data has been normalised using the friction velocity, u τ , and the mean channel half-height, δ. A slight transparency has been applied to the contours to show the underlying roughness topography.

T.O. Jelly and A. Busse
International Journal of Heat and Fluid Flow 80 (2019) 108485 occupied" region become comparable, i.e. u ũ 1 1 . Moltchanov and Shavit (2013) combined this reasoning with a parameterisation of the wake-occupied area and devised a phenomenological closure model to predict D 11 . They tested the performance of their model using high-resolution PIV data acquired inside and around randomly distributed thin vertical glass plates and cylinders. Overall, a good level of agreement between the modelled and measured stresses was reported, with mean absolute errors quoted on the order of 10%. In what follows below, the Moltchanov and Shavit (2013) model will be tested on the current data to establish whether similar levels of agreement can be obtained for an irregular Gaussian surface. Moltchanov and Shavit (2013) proposed that the streamwise dispersive stress, D 11 , and the square of the DA streamwise velocity, u , 1 2 are linearly related through a model parameter they called the "relative wake area". The relative wake area can be expressed as where A w is the wake-occupied area on each wall-parallel plane. Following Moltchanov and Shavit (2013), wake-occupied regions are defined as those where the time-averaged streamwise velocity is less than half the local DA value, i.e. fluid-occupied points that satisfy the con- 3 . Before proceeding, it is worth noting that the concept of a "wake-occupied region" needs to be interpreted with care in the context of the irregular surface under consideration in this work. This is because some of the flow below the mean roughness plane (x 3 /δ < 0) recirculates within the roughness cavities (Fig. 8), and bears no obvious resemblance to a wake. In contrast, wake-like features are clearly visible in the upper half of the roughness canopy where localised patches of low-momentum flow occur immediately downstream of the protruding peaks (see Fig. 13(a)-(e)). Therefore, we anticipate the model devised by Moltchanov and Shavit (2013) to perform best at wall-normal heights that coincide with these "wake-occupied" regions. Moltchanov and Shavit (2013) proposed two different models using the relative wake area (Eq. (15)). The first model assumes that the level of dispersive stress is directly proportional to the local value of the relative wake area, which can be expressed as where subscript "m1" denotes the first of the two models. Note that model m1 only takes into account wake-occupied regions, where the inequality 3 is satisfied. In their second model, Moltchanov and Shavit (2013) split the timeaveraged velocity field into two separate parts: (i) The wake-occupied part, which they defined as the α fraction of the total domain area and (ii) the wake-unoccupied part defined by the remaining (1 ) fraction. By evaluating the streamwise dispersive velocity in regions (i) and (ii), Moltchanov and Shavit (2013) obtained the following expression where subscript "m2" denotes the second of the two dispersive stress models. Before comparing the modelled and simulated stresses, the constitutive ingredients of model m1 (Eq. (16)) and model m2 (Eq. (17)) will be examined in further detail. The first ingredient of the Moltchanov and Shavit (2013) model is the relative wake area (Eq. (15)). Wall-normal profiles of the relative wake area, α, are plotted for each friction Reynolds number in Fig. 14. As the highest roughness crest is approached, the relative wake area decays towards zero -implying that negligible levels of dispersive stress will be predicted in this region and in the outer flow. Below the crests, the relative wake area begins to increase and varies linearly on the range 0 ≲ x 3 /h max ≲ 0.5 before attaining its maximum value of α ≈ 0.6 about the roughness mean plane. The inset panel on Fig. 14 shows profiles of /(1 ) for each friction Reynolds number. Whilst the profiles of /(1 ) show the same Reynolds-number trends as the profiles of α, their magnitude is approximately a factor of two higherwhich means that the second model (Eq. (17)) will predict stronger dispersive stresses than that of the first (Eq. (16)).
The second ingredient of the Moltchanov and Shavit (2013) model is the square of the DA streamwise velocity, + u 1 2 . Profiles of the squared velocity profiles for each friction Reynolds number are plotted in Fig. 15. Above the highest roughness crest, the profiles of + u 1 2 decrease with increasing friction Reynolds number -consistent with the inner-and outer-scaled velocity profiles shown previously in Figs. 6 and 7. Beneath the roughness crests, the magnitude of + u 1 2 increases with increasing friction Reynolds number on the range 0 ≲ x 3 /h max ≲ 0.6. Considering the relative magnitude and Reynolds number insensitivity of the relative wake area in the same region (see Fig. 14), an increase in + u 1 2 will lead to higher levels of dispersive stress. The inset panel of Fig. 15 shows a zoomed-in view of the bottom half of the roughness canopy. As expected, the profiles of + u 1 2 in this region follow the same Reynolds number trends as the mean reverse profiles shown previously in Fig. 7.   Fig. 14. Profiles of relative wake area, α. Line-types are specified in Table 3. All data has been normalised using the friction velocity, u τ , and the height of the highest roughness crest, h max . The height of the highest roughness crest is shown as a vertical gray line. The inset panel shows wall-normal profiles of the ratio /(1 ).

Fig. 15.
Profiles of the squared DA streamwise velocity, + u , 1 2 in the nearroughness region. Line-types are specified in Table 3. All data has been normalised using the friction velocity, u τ , and the height of the highest roughness crest, h max . The height of the highest roughness crest is shown as a vertical gray line. The inset panel shows a zoomed in view of the squared DA velocity profiles in the bottom half of the roughness canopy.
T.O. Jelly and A. Busse International Journal of Heat and Fluid Flow 80 (2019) 108485 The performance of the dispersive stress models was quantified by computing the mean absolute error (MAE), defined here as where subscript "mi" denotes a quantity modelled using either model m1 (Eq. (16)) or model m2 (Eq. (17)). The MAE and PSE on the modelled stresses for each friction Reynolds number are compared in Table 6. The error analysis indicates that model m2 consistently outperforms model m1 for the friction Reynolds numbers range of this study -which agrees with the past analysis of Moltchanov and Shavit (2013). Relative to the target DNS data, the MAE corresponding to model m2 is, at the very worst, on the order of 10%. In addition, the PSE corresponding to model m2 matches the DNS data to within approximately 5% for all friction Reynolds numbers considered in this study.
Profiles of the streamwise dispersive stress predicted using the second model (Eq. (17)) are compared against the DNS data in Fig. 16. Overall, the shape and magnitude of the predicted and simulated stress profiles are in good agreement for each value of the friction Reynolds number. In particular, the peak magnitude of dispersive stress and its wall-normal position are both captured with an impressive level of accuracy. In addition, negligible levels of dispersive stress are also correctly predicted in the region above the roughness crests (x 3 / h max > 1). In contrast, the model struggles to reproduce the levels of dispersive stress in the vicinity of the canopy upper edge (x 3 /h max ≈ 1) and in the lower half of the roughness canopy The mismatch between the modelled and simulated stresses in the upper and lower regions of the roughness canopy highlight some key differences between this work and the past study of Moltchanov and Shavit (2013). First, as was previously mentioned, the concept of a "wake-occupied region" in the lower half of the roughness canopy cannot be applied to the current surface. As a result, neither model m1 (Eq. (16)) or model m2 (Eq. (17)) can be expected to perform well in this region. Second, whereas the void fraction under consideration in this study varies smoothly throughout the roughness canopy (see Fig. 5), the same quantity remains constant beneath the crests of the constant cross-section canopies studied by Moltchanov and Shavit (2013). Third, whereas the void fraction of the irregular surface varies smoothly across the roughness crests, the same quantity becomes discontinuous for the idealised vegetation canopies studied by Moltchanov and Shavit (2013). Sudden jumps in the void fraction function will lead to sharp changes in the relative wake area (Eq. (15)), ultimately leading to stronger levels of streamwise dispersive stress in upper region of the canopy. Such behaviour is clearly evident in the past results of Moltchanov and Shavit (2013) (e.g. see Figs. 5-7 in their paper). A different behaviour is observed in this study, whereby the void fraction function (Fig. 5), the relative wake area (Fig. 14) and the dispersive stresses (Fig. 16) all exhibit slowly increasing and smoothly varying profiles just below the roughness crests, before reaching their maxima in the upper half of the canopy. Despite these differences, the agreement between the modelled and simulated stress profiles is still quite good (see Fig. 16 and Table 6). Therefore, although Moltchanov and Shavit (2013) formulated their model (Eqs. (16) and (17)) specifically in the context of idealised vegetation canopies, the results of this study imply that the framework they developed is sufficiently robust to predict D 11 above and within the wake-occupied regions of an irregular near-Gaussian surface.
Finally, it is worth noting that the prediction of D 11 using Eqs. (16) and (17) requires prior knowledge of the DA velocity profile, u 1 . Moreover, evaluating the relative wake area (Eq. (15)) requires prior knowledge of the time-averaged streamwise velocity field, u , 1 including in the region below the roughness crests. Whilst such information can be obtained directly from DNS data (see Fig. 13), acquiring an equivalent dataset in the laboratory is very challenging -especially for irregular roughness topographies where optical access may be limited. Ideally, an estimate for the relative wake area could be obtained a priori for an arbitrary roughness geometry. Some progress towards this goal was made by Moltchanov and Shavit (2013), who calculated the relative wake area of their cylinder canopy to within approximately 15% of the measured value based on past empirical relationships devised by Zdravkovich (1997). However, as noted by the authors, their estimate of the relative wake area was only valid provided that the wake width was equal to the cylinder diameter and that wake-wake interactions remained negligible. With reference to Fig. 13, extrapolating such analysis to the irregular roughness topography under consideration in this work is clearly not straightforward, and, as a result, more data for a variety of surface topographies and flow conditions is needed.

Summary and conclusions
DNS of turbulent channel flow with irregular near-Gaussian roughness on both the upper and lower walls (Fig. 3) were performed at four separate friction Reynolds numbers, namely, 180, 240, 360 and 540. In order of increasing Re τ , the roughness functions obtained from each simulation were found to be = + U {4.27, 5.33, 6.75, 8.17}. Based on an inferred value of Nikuradse's equivalent sandgrain roughness, the data included in this study spans the transitionally and fully rough regimes (Fig. 4). The Reynolds-number dependence of DA statistics was studied in the context of streamwise velocity profiles (Figs. 6 and 7), mean recirculation regions (Fig. 8), dispersive and Reynolds stresses  Table 3. Symbols correspond modelled stress data at = × Re {180( ), 240( ), 360( ), 540( )}. All data has been normalised using the friction velocity, u τ , and the height of the highest roughness crest, h max . The height of the highest roughness crest is shown as a vertical gray line.

T.O. Jelly and A. Busse
International Journal of Heat and Fluid Flow 80 (2019) 108485 (Figs. 9 and 10), an analysis of the mean force balance above the roughness crests (Fig. 11) and a comparison of the shear stress gradients (Fig. 12), as well as time-averaged stress distributions within the roughness canopy (Fig. 13). In addition, the constitutive ingredients of the Moltchanov and Shavit (2013) model (Eqs. (16) and (17)) were examined (Figs. 14 and 15), and, finally, a comparison between the modelled and simulated streamwise dispersive stresses was made (Fig. 16).
The principal aim of this study was to study the dependence of DA velocity and stresses as a function of the friction Reynolds number for an irregular Gaussian surface. This work was motivated by several recent studies that adopted the opposite approach, i.e. those that considered multiple roughness geometries at a fixed friction Reynolds number. The key findings of this study are summarised below.
(i) For the friction Reynolds number range considered in this study, the DA velocity and Reynolds stress profiles are in strong support of Townsend's outer-layer similarity hypothesis (Townsend, 1976), e.g. see Figs. 6, 9 and 10. However, despite the collapse of the smooth-and rough-wall data in the outer region, this work highlights the difficulty in defining the length-scale ratio δ/k for surfaces with a broad range of topographical scales. For instance, well-established criteria, such as the threshold of δ/k > 40 for the preservation of wall similarity above transverse square bar roughness (Jiménez, 2004), become ambiguous for irregular surfaces. This is because multiple measures of k exist for a single heightmap (see Table 4). In addition, for each measure of k there is a corresponding value of + k , which, in turn, complicates the classification of roughness effects based on the "traditional" sandgrain approach. For the irregular Gaussian surface in this work, the mean peak-to-valley was found to match the equivalent sandgrain roughness to within 10% (see Fig. 4). This implies that S z,5 × 5 is an appropriate hydraulic length-scale to predict roughness effects in the fully-rough regime (for this particular roughness) -supporting the past work by Busse et al. (2017) and Thakkar et al. (2018). (ii) Above the roughness crests, the dispersive stresses become negligible for each friction Reynolds number considered in this study (see Figs. 9 and 10). This finding agrees with recent work regarding turbulent flow past irregular, isotropic distributions of small-scale roughness (Forooghi et al., 2018) but differs from other studies that consider surfaces with highly directional, larger-scale features. For instance, past experimental work by Barros and Christensen (2014) reports secondary rolls in the time-averaged velocity field above a streamwise anisotropic irregular surface, that, upon double-averaging, would yield non-negligible dispersive stresses in the outer region. Beneath the crests of the irregular Gaussian surface, the levels of dispersive stress grow with increasing Reynolds number, and, for the highest value of Re τ considered in this work, each component of "form-induced" stress exceeds its Reynolds stress counterpart. The streamwise component of dispersive stress is always dominant, followed by the spanwise and wall-normal stresses -indicating that the strongest deviations from streamwise-spanwise homogeneity occur along the primary flow direction. (iii) In the context of mean flow dynamics, the relative importance of dispersive and Reynolds shear stresses is determined by the respective magnitude (and sign) of their gradient. This is because both terms appear on the right-hand side of the mean force balance Eq. (14) and their combined effect maintains the flow against the mean pressure gradient and viscous force above the roughness crests (see Fig. 11). Just above the crests, the dispersive and Reynolds shear stress gradients make an approximately equal (and opposite) contribution to the mean force balance, where the latter serves as a momentum sink and and former serves as a momentum source. A similar source / sink scenario has been described in the context of fully-developed turbulent flow past idealised vegetation canopies (Moltchanov et al., 2015). In contrast, comparing the shear stress magnitudes, as opposed to their gradient, implies that the level of dispersive shear stress is negligible at the roughness crests, relative to Reynolds shear stress. Therefore, in-line with the past work of Manes et al. (2008), our analysis shows that it might be misleading to quantify the dynamical significance of either shear stress using just its magnitude. (iv) Dispersive stresses represent the mean-squared signature of stationary flow patterns, which, in this study, appear predominantly in the form of "wake-occupied regions" downstream of the roughness peaks (see Fig. 13). Strong streamwise dispersive fluxes develop in these regions because the difference between the timeaveraged and DA averaged velocity is large. The opposite behaviour is observed in the spaces between the peaks, i.e. the "wakeunoccupied" region, where the DA and time-averaged velocities are approximately equal. Furthermore, regions of strong streamwise dispersive stress correlate with weak streamwise Reynolds stress and vice versa. As the friction Reynolds number increases, the levels of dispersive stress in wake-occupied regions also increases, which, ultimately, leads to larger DA values (see Fig. 9(a) and (b)). (v) The levels of streamwise dispersive stress were predicted to within 10% of the DNS results using a model devised by Moltchanov and Shavit (2013). This model assumes that levels of dispersive stress scale with the square of the DA velocity, i.e. D u , 11 1 2 and relies on a parameterisation of the wake-occupied area -the so-called "relative wake area" (Eq. (15)) -to predict the levels of "forminduced" stress above (and within) the roughness canopy. The suitability of this model beneath the mean roughness plane is questionable because the flow in this region and bears little resemblance to a wake. Nonetheless, good levels of agreement above the roughness mean plane were observed, and, in particular, the peak value of dispersive stress (and its position) was captured with a remarkable level of accuracy for each value of the friction Reynolds number considered in this study (see Fig. 16). (vi) Despite a good match between the simulated and modelled streamwise dispersive stresses (Fig. 16), it is important to recognise that term D 11 makes zero contribution to the momentum balance in a fully-developed rough-wall turbulent channel flow (Eq. (14)). In contrast, the streamwise gradient of D 11 can play an important role in determining the mean force balance of spatially-developing flows, e.g. see the order-of-magnitude analysis in past work by Moltchanov et al. (2011). Related to this point is how to predict the remaining dispersive normal stresses, and, more importantly, the dispersive shear stress, D 13 , whose wall-normal gradient plays a key role in maintaining the balance of mean forces above the roughness crests (Fig. 11). However, as was previously mentioned, there does not appear to be an obvious connection between the "wake-occupied" regions and the generation of dispersive shear stress (Fig. 13) -implying that an accurate prediction of D 13 cannot be obtained via the original Moltchanov and Shavit (2013) model. This observation is supported by the past work of Finnigan et al. (2015), who noted that adopting a "velocitysquared closure" (such as Eq. (16) and (17)) can fail to predict the levels of dispersive shear stress even for simple canopy flows. As a result, more traditional approaches, such as those based on eddy diffusivity and mixing length arguments, may give better results.
Finally, whereas this study has focused on the fluid dynamic properties of an irregular near-Gaussian surface, naturally-occurring surfaces often exhibit highly non-Gaussian properties, e.g. see the skewness-kurtosis map plotted in Fig. 1. A number of past experimental (Schultz and Flack, 2009;Flack and Schultz, 2010) and numerical (Forooghi et al., 2018;Jelly and Busse, 2018b) studies have identified skewness as a key topographical parameter in the context of rough-wall T.O. Jelly and A. Busse International Journal of Heat and Fluid Flow 80 (2019) 108485 turbulent flows. As a result, future work should focus on the Reynolds number dependence of surfaces with systematically varied level of skewness with the ultimate aim of predicting roughness effects in application-scale flow configurations. To this end, performing an extended literature survey to determine the typical skewness-kurtosis combinations encountered in practice may prove an useful first step (see Fig. 1).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.