1. Introduction
Turbulent flows bounded by rough walls are abundant in both nature – e.g. fluvial flows (Mazzuoli & Uhlmann Reference Mazzuoli and Uhlmann2017) and wind flow over vegetation (Finnigan & Shaw Reference Finnigan and Shaw2008) and urban canopies (Coceal & Belcher Reference Coceal and Belcher2004; Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016) – and industry – e.g. degraded gas turbine blades (Bons et al. Reference Bons, Taylor, McClain and Rivir2001), bio-fouled ship hulls (Hutchins et al. Reference Hutchins, Monty, Nugroho, Ganapathisubramani and Utama2016), iced surfaces in aero-engines (Velandia & Bansmer Reference Velandia and Bansmer2019) and deposited surfaces inside combustion chambers (Forooghi et al. Reference Forooghi, Weidenlener, Magagnato, Böhm, Kubach, Koch and Frohnapfel2018c). Systematic study of roughness effects on skin friction dates back to the pioneering works of Nikuradse (Reference Nikuradse1933) and Schlichting (Reference Schlichting1936). Flow-related roughness is usually classified into two types,
$k$-type and
$d$-type roughness. For
$k$-type roughness, the flow response depends directly on the physical scale of the roughness height; while for
$d$-type roughness, it is insensitive to the roughness height scale, but rather determined by the outer length scales, e.g. the pipe diameter
$d$ (Perry, Schofield & Joubert Reference Perry, Schofield and Joubert1969; Jiménez Reference Jiménez2004). The scope of the present work is limited to three-dimensional (3-D) irregular rough surfaces, in which
$k$-type behaviour is relevant.
In industry, the Moody diagram (Moody Reference Moody1944) has been considered as a standard method to calculate the skin friction of a rough surface. The Moody diagram relates the friction factor to the roughness height
$\epsilon$, which is linked to the equivalent sand-grain roughness
$k_s$. This quantity, which is also an input to many low-fidelity turbulence models for rough walls (Suga, Craft & Iacovides Reference Suga, Craft and Iacovides2006; Brereton & Yuan Reference Brereton and Yuan2018), is not known a priori for any given irregular rough surface. Hence, for any new roughness topography,
$k_s$ needs to be determined using laboratory or high-fidelity numerical experiments or estimated based on roughness correlations derived from such experiments.
The problem of predicting the roughness-induced friction drag based merely on knowledge of the roughness topography has received extensive attention in the past, and a variety of roughness correlations have been developed in different industrial contexts (Waigh & Kind Reference Waigh and Kind1998; Macdonald Reference Macdonald2000; van Rij, Belnap & Ligrani Reference van Rij, Belnap and Ligrani2002; Bons Reference Bons2005; Flack & Schultz Reference Flack and Schultz2010; Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2015; Forooghi et al. Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017; Thakkar, Busse & Sandham Reference Thakkar, Busse and Sandham2017; Flack, Schultz & Barros Reference Flack, Schultz and Barros2020). In these roughness correlations, the topography of the rough surface is often represented by statistical measures of the roughness height map
$k(x,z)$, with
$k$ being the surface height as a function of horizontal coordinates
$x$ and
$z$. Some widely discussed statistical measures in this context are summarized in the recent review of the topic by Chung et al. (Reference Chung, Hutchins, Schultz and Flack2021), for instance the skewness
$Sk$ (Flack & Schultz Reference Flack and Schultz2010; Forooghi et al. Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017), effective slope
$ES$ (Napoli, Armenio & De Marchis Reference Napoli, Armenio and De Marchis2008; Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2015) and density parameter
$\varLambda _s$ (Sigal & Danberg Reference Sigal and Danberg1990; van Rij et al. Reference van Rij, Belnap and Ligrani2002). Despite extensive work in the past, a universal correlation with the ability to accurately predict the drag of a generic rough surface remains elusive (Flack Reference Flack2018). Arguably, the development of such a correlation requires a large amount of data from realistic roughness samples. However, the generation of an appropriate database has been hindered mainly due to two factors: the formidable cost associated with many numerical or laboratory experiments, and the relative scarceness of realistic roughness maps combined with the lack of ability to systematically vary their properties.
A considerable portion of data in the literature deals with regular roughness – often generated by the distribution of similar geometric elements. Examples of the geometries studied include cubes (Orlandi & Leonardi Reference Orlandi and Leonardi2006; Leonardi & Castro Reference Leonardi and Castro2010), spheres (Mazzuoli & Uhlmann Reference Mazzuoli and Uhlmann2017), pyramids (Schultz & Flack Reference Schultz and Flack2009), LEGO bricks (Placidi & Ganapathisubramani Reference Placidi and Ganapathisubramani2015), ellipsoidal egg-carton shape (Bhaganagar Reference Bhaganagar2008) and sinusoidal roughness (Chan et al. Reference Chan, MacDonald, Chung, Hutchins and Ooi2015, Reference Chan, MacDonald, Chung, Hutchins and Ooi2018). In comparison, investigations based on realistic rough surfaces are less frequent and include a much lower number of cases. Notably, Thakkar et al. (Reference Thakkar, Busse and Sandham2017) utilized direct numerical simulation (DNS) to study the effect of roughness topography on flow statistics for 17 industrially relevant irregular surfaces, and proposed roughness correlations for the transitionally rough regime. Other examples of realistic roughness studies in the framework of wall-bounded turbulence include Cardillo et al. (Reference Cardillo, Chen, Araya, Newman, Jansen and Castillo2013), Yuan & Piomelli (Reference Yuan and Piomelli2014), Busse, Lützner & Sandham (Reference Busse, Lützner and Sandham2015), Busse, Thakkar & Sandham (Reference Busse, Thakkar and Sandham2017), Forooghi et al. (Reference Forooghi, Weidenlener, Magagnato, Böhm, Kubach, Koch and Frohnapfel2018c), Yuan & Jouybari (Reference Yuan and Jouybari2018), Jouybari, Brereton & Yuan (Reference Jouybari, Brereton and Yuan2019) and Mangavelli, Yuan & Brereton (Reference Mangavelli, Yuan and Brereton2021).
In recent years, mathematically generated surfaces have received an increased amount of attention as a means to systematically study realistic irregular roughness. Many of these roughness generation approaches rely on random superposition of discrete geometric elements (Scotti Reference Scotti2006; Chau & Bhaganagar Reference Chau and Bhaganagar2012; Forooghi, Stripf & Frohnapfel Reference Forooghi, Stripf and Frohnapfel2018a; Kuwata & Kawaguchi Reference Kuwata and Kawaguchi2019) or Fourier modes (Anderson & Meneveau Reference Anderson and Meneveau2011; De Marchis, Saccone & Milici Reference De Marchis, Saccone and Milici2020). Some authors (Barros, Schultz & Flack Reference Barros, Schultz and Flack2018; Jelly & Busse Reference Jelly and Busse2019) opt for approaches based on the linear combination of random numbers in an attempt to study surfaces that resemble realistic roughness as closely as possible. These methods can produce a prescribed power spectrum (PS) with a Gaussian probability density function (p.d.f.) for the generated roughness. Motivated by the relevance of non-Gaussian roughness in industry, recently Flack et al. (Reference Flack, Schultz and Barros2020) employed a modified version of such methods to study non-Gaussian roughness with a certain choice of p.d.f.s. In the present paper, we adopt an alternative roughness generation method proposed by Pérez-Ràfols & Almqvist (Reference Pérez-Ràfols and Almqvist2019). This method is deemed advantageous in producing surrogates of realistic roughness since it provides absolute flexibility to prescribe any desired combination of p.d.f. and PS, as well as more robustness compared to algorithms based on translations of the Pearson's or Johnson's types used by previous authors. We refer to the roughness samples generated by this method as ‘pseudo-random’ roughness, in the sense that the topography is random but its statistical properties are prescribed.
In recent years, DNS has been the pacing approach in studying the effect of roughness topography on friction drag. Standard DNS, however, involves resolving the entire spectrum of turbulent length scales ranging from large geometrical scales to the small viscous scale, which is computationally costly. To tackle this problem, Chung et al. (Reference Chung, Chan, MacDonald, Hutchins and Ooi2015) and MacDonald et al. (Reference MacDonald, Chung, Hutchins, Chan, Ooi and Garcıa-Mayoral2016) employed the idea of DNS in minimal span channels (Jiménez & Moin Reference Jiménez and Moin1991; Flores & Jiménez Reference Flores and Jiménez2010) for prediction of roughness-induced drag over a regular sinusoidal roughness in a channel. The central idea followed by these authors is that the amount of downward shift in the inner-scaled velocity profile
${\rm \Delta} U^+$ is the determining factor in the prediction of drag. These authors showed that, thanks to outer-layer similarity of wall-bounded turbulence (Townsend Reference Townsend1976), this key quantity can be accurately predicted by minimal rough channels. This can be achieved as long as an adequately large range of near-wall structures are accommodated in the simulation domain. The same group of authors further developed the idea for channels with minimal streamwise extent (MacDonald et al. Reference MacDonald, Chung, Hutchins, Chan, Ooi and Garcia-Mayoral2017), for high-aspect-ratio transverse bars (MacDonald et al. Reference MacDonald, Ooi, García-Mayoral, Hutchins and Chung2018) and also for passive scalar calculations (MacDonald, Hutchins & Chung Reference MacDonald, Hutchins and Chung2019).
These efforts established the following criteria for the size of a minimal channel based on simulations with 3-D sinusoidal roughness:
\begin{equation} L_z^+\geq \text{max}\left(100,\frac{\tilde{k}^+}{0.4},\lambda_{sin}^+\right),\quad L_x^+\geq \text{max}(1000,3L_z^+,\lambda_{sin}^+). \end{equation}
Here
$L_z$ and
$L_x$ are the spanwise and streamwise extents of the minimal channel, respectively;
$\lambda _{sin}$ is the sinusoid wavelength of roughness;
$\tilde {k}$ is the characteristic roughness height (here sinusoid amplitude); and the
$+$ superscript indicates viscous scaling. While the aforementioned studies showed the potential of minimal channels in the determination of roughness-induced drag, a formal extension of this concept to random, irregular roughness is yet to be made. Recently, Alves Portela & Sandham (Reference Alves Portela and Sandham2020) reported predictions of flow over a realistic roughness, combining minimal channels with a novel hybrid direct numerical simulation / unsteady Reynolds-averaged Navier–Stokes (DNS/URANS) model. These authors highlighted the need for careful investigation of the minimal-channel concept for realistic roughness. Although some previous researchers, e.g. Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021) and Pargal, Yuan & Brereton (Reference Pargal, Yuan and Brereton2021), have already applied the concept successfully for random surfaces, no systematic verification of the minimal-channel approach for irregular roughness of random nature has been reported in the past.
In view of the above, the present work aims to provide a systematic proof for the validity of the minimal-channel approach in irregular roughness of random nature. Here ‘systematic’ refers to covering a wide range of roughness topographies (12 different p.d.f./PS combinations) and both transitionally and fully rough regimes. The study aims to answer two key questions. The first is whether it is adequate for roughness in the minimal channel to be similar to the original roughness merely in a ‘statistical’ sense. This, in turn, relates to the more fundamental question about the minimum amount of statistical information needed to uniquely predict the roughness-induced drag. The second question is whether the rules set by (1.1a,b) can be modified or relaxed so that they are applicable to any type of roughness. Here the critical issue is that a realistic roughness may contain very large horizontal wavelengths, making the minimal-channel approach futile if the original rules are to be met strictly. We compare the results from minimal- and full-channel DNS in § 3.1. To shed more light on the latter question, in § 3.2 we also study in detail the local distribution of drag force on a rough surface to better understand the contribution of different horizontal scales in drag generation. As a final point, we use the generated data to assess a number of widely used roughness correlations in § 3.3. A detailed description of our methodology and a summary of findings are presented in §§ 2 and 4, respectively.
2. Numerical methodology
2.1. Pseudo-random roughness generation
As mentioned in the introduction, the roughness generation method proposed by Pérez-Ràfols & Almqvist (Reference Pérez-Ràfols and Almqvist2019) is adopted in the present methodology. In this method both the wall-parallel and the wall-normal statistical properties of the roughness can be adjusted. Here wall-parallel properties refer to the PS of the roughness structures, and wall-normal properties refer to the p.d.f. of the surface height. The roughness map is represented by a discrete elevation distribution on a two-dimensional (2-D) Cartesian grid. The generation algorithm used in the present work takes the target p.d.f. and PS as inputs. Transformations between the physical space and spectral space are done by a discrete fast Fourier transform. Initially a roughness map
$k_{pdf}^0$ is generated that has the prescribed p.d.f. but not necessarily the prescribed PS. This initial map is then corrected in the Fourier space according to the prescribed PS, which is represented by
$\hat {k}_{{PS}}^{i}$:
\begin{equation} \hat{k}_{{PS}}^{i+1} =\hat{k}_{pdf}^{i}\frac{|\hat{k}_{{PS}}^{i}|}{|\hat{k}_{pdf}^{i}|}, \end{equation}
where
$i$ indicates the iteration of the generation process. The output of this stage,
$k_{{PS}}^{i+1}$, has the desired PS but not necessarily the prescribed p.d.f. In the present notation, subscripts
$pdf$ and
$PS$ indicate that the roughness field has the desired p.d.f. or PS, respectively. The hat indicates the Fourier transform. After which,
$k_{pdf}^{i+1}$ is updated by correcting the p.d.f. of
$k_{{PS}}^{i+1}$ by rank ordering. This correction process continues for
$n$ iterations until both
$k_{pdf}^n$ and
$k_{{PS}}^n$ converge to a height map with the target p.d.f. and PS within a predetermined error. For more details on the generation algorithm, readers are referred to Pérez-Ràfols & Almqvist (Reference Pérez-Ràfols and Almqvist2019).
2.2. Direct numerical simulation
A number of DNS have been carried out in fully developed turbulent plane channels, in which the flow is driven by a constant pressure gradient. A representation of the simulation domain is shown in figure 1, where
$x$,
$y$ and
$z$ denote the streamwise, wall-normal and spanwise directions, with respective velocity components
$u$,
$v$ and
$w$. The roughness elements are mounted on both the upper wall and the lower wall. The channel half-height, which is the distance between the deepest trough in the roughness and the centre-plane of the channel, is labelled as
$H$, and will be used to normalize the geometrical scales. The incompressible Navier–Stokes equations are solved using the pseudo-spectral solver SIMSON (Chevalier et al. Reference Chevalier, Schlatter, Lundbladh and Henningson2007), where wall-parallel directions are discretized in Fourier space, while in the wall-normal direction Chebyshev discretization is employed. The immersed boundary method (IBM) based on Goldstein, Handler & Sirovich (Reference Goldstein, Handler and Sirovich1993) is used to impose the no-slip boundary condition on the roughness by introducing an external volume force field directly to the Navier–Stokes equation. The presently used IBM implementation has been validated and used in previous studies (Forooghi et al. Reference Forooghi, Stroh, Schlatter and Frohnapfel2018b; Vanderwel et al. Reference Vanderwel, Stroh, Kriegseis, Frohnapfel and Ganapathisubramani2019; Stroh et al. Reference Stroh, Schäfer, Frohnapfel and Forooghi2020).
Schematic representation of the simulation domain with an example pseudo-realistic surface mounted. Normalization of lengths with
$H$ is applied in this figure.

The Navier–Stokes equation can be written as
where
$\boldsymbol {u}=(u,v,w)^{\rm T}$ is the velocity vector and
$P_x$ is the mean pressure gradient in the flow direction added as a constant and uniform source term to the momentum equation to drive the flow in the channel. Here
$p$,
${\boldsymbol {e}}_x$,
$\rho$,
$\nu$ and
$\boldsymbol {f}_{IBM}$ are pressure fluctuation, streamwise basis vector, density, kinematic viscosity and external body-force term due to IBM, respectively. Periodic boundary conditions are applied in the streamwise and spanwise directions. A no-slip boundary condition is applied on the rough walls. The friction Reynolds number is defined as
$Re_\tau =u_\tau (H-k_{{md}})/\nu$, where
$u_\tau =\sqrt {\tau _w/\rho }$ and
$\tau _w=-P_x (H-k_{{md}})$ are the friction velocity and wall shear stress, respectively. The meltdown height, denoted by
$k_{{md}}$, is the mean roughness height measured from the deepest trough. Note that we use
$(H-k_{{md}})$, which is the mean half-cross-sectional area divided by the channel width – or arguably the effective channel half-height – as the reference length scale in these definitions. In total, four different values of
$Re_\tau$ in the range of 250–1000 are simulated in the present work in order to be able to cover different regimes. The simulations designed to study the effect of roughness topography are, however, performed mainly at
$Re_\tau =500$.
The simulation domain is discretized in an equidistant grid in wall-parallel directions, while in the wall-normal direction cosine stretching based on Chebyshev node distribution is applied. The selection of grid size must take into consideration both flow and roughness length scales. As reflected by Busse et al. (Reference Busse, Lützner and Sandham2015), each roughness wavelength should be represented by multiple computational cells. Since we prescribe the PS in the roughness generation approach, the range of present roughness wavelengths can be prescribed. Here the smallest roughness wavelength, labelled as
$\lambda _1$, is the crucial quantity for horizontal grid resolution. Therefore, in view of the present computational capacity,
$\lambda _1=0.08H$ is prescribed for the following simulations. A mesh independence study is carried out, from which the smallest roughness wavelength being resolved by 8–10 cells in each direction is found to be adequate to obtain the converged double-averaged velocity profile. Overall, the grid size
$\varDelta ^+<5$ in wall-parallel directions is proven to be appropriate through the mesh independence test. In the wall-normal direction, a cosine stretching mesh is adopted for the Chebyshev discretization. It is also checked through the mesh independence test that, for the present types of rough surfaces, a vertical cell number of 401 is sufficient to deliver a converged result at the highest
$Re_\tau \approx 1000$, thanks to the overresolving of the roughness structure by the cosine stretching grid near the wall.
2.3. Description of cases
Using the roughness generation algorithm introduced in § 2.1, multiple samples are generated. In the present work, different types of p.d.f. will be combined with power-law PS, that is,
$E_k(\boldsymbol {q})=C_0(\|\boldsymbol {q}\|/q_0)^p$, where
$\boldsymbol {q}=(q_x,q_z)^{\rm T}$ is the wavenumber vector,
$q_0=2{\rm \pi} /\lambda _0$ is the smallest wavenumber corresponding to the largest in-plane length scale
$\lambda _0$,
$C_0$ is a constant to scale the roughness height, and
$p$ is the spectral slope of the power-law PS. An overview of the configurations of p.d.f. and PS is illustrated in figure 2. In figure 2(b–d) the PS density normalized by the root mean square (r.m.s.) of the roughness height are compared in pairs. The upper and lower cutoff wavelengths
$\lambda _0$ and
$\lambda _1$ are transformed to cutoff wavenumbers
$q_0=2{\rm \pi} /\lambda _0$ and
$q_1=2{\rm \pi} /\lambda _1$, which are represented by the red dashed lines in figure 2(b–d) on the left and right sides, respectively. As stated in the previous section, the lower cutoff wavelength is related to the grid resolution, and a value of
$\lambda _1=0.08H\approx 8\varDelta _{x}\approx 8\varDelta _{z}$ is applied for all roughness topographies in the present work. With an isotropic roughness and a fixed
$\lambda _1$, the PS is determined by two remaining parameters,
$\lambda _0$ and
$p$. In the present work, two values of
$p$ (
$p=-1$ and
$p=-2$) are examined, the PS of which are shown in figure 2(b). For the selection of
$p$ values, we seek similarity to previous works (Anderson & Meneveau Reference Anderson and Meneveau2011; Barros et al. Reference Barros, Schultz and Flack2018; Nikora et al. Reference Nikora2019). Moreover, two different upper cutoff wavelengths (
$\lambda _0=0.8H$ and
$\lambda _0=1.6H$) of the roughness PS are investigated. The PS with
$\lambda _0=0.8H$ and
$1.6H$ with identical slopes
$p$ are compared in figure 2(c,d), where wavenumber
$q$ is normalized by referencing wavenumber
$q_{{ref}}=2{\rm \pi} /(0.8H)$
Statistical representation of the studied roughness. (a) The p.d.f. of roughness, (b) normalized PS density with different
$p$, (c) normalized PS density with different
$\lambda _0$, for
$p=-1$, and (d) normalized PS density with different
$\lambda _0$, for
$p=-2$. In (b–d) wavenumber
$q$ is normalized by the reference wavenumber
$q_{{ref}}=2{\rm \pi} /(0.8H)$. Vertical dashed lines are high-pass filtering and low-pass filtering, corresponding to
$\lambda _0$ and
$\lambda _1$, respectively.

Three types of p.d.f.s with positive, zero and negative skewness are examined in the present work. Covering a relatively large range of
$Sk$ is intended to ensure that the results can be generalized to a wide spectrum of naturally occurring roughness in different applications. The non-skewed roughness is described by a Gaussian distribution. For the positively skewed roughness, a Weibull distribution is used, which can be written as
where the parameters
$K$ and
$\beta$ can be used to adjust the standard deviation and skewness of the distribution. The skewness is always adjusted to the value of 0.48. Similar to the Gaussian distribution, the kurtosis of the Weibull distribution is always equal to 3. A negatively skewed p.d.f. is obtained by flipping the p.d.f. of a positively skewed Weibull p.d.f. Here
$Sk=-0.48$ is prescribed. In the present work, the
$99\,\%$ confidence interval of roughness height p.d.f.,
$k_{99}$, is used as the characteristic size of roughness, i.e.
$k=k_{99}$. This measure is related to the standard deviation of the roughness, and hence can be directly prescribed. We used a fixed value of
$k=0.1H$ in all cases. Indeed,
$k$ is a statistical measure of the maximum peak-to-trough roughness size, which, unlike the absolute peak-to-trough size
$k_{t}$, is not deteriorated by extreme events. These types of p.d.f.s are illustrated in figure 2. Moreover, in order to avoid extreme high roughness elevations in the simulations, roughness heights outside 1.2 times the
$99\,\%$ confidence interval of the p.d.f. are excluded. Therefore, the peak-to-trough height
$k_{t}=0.12H$ is achieved for all roughness in the present work. Combining the three types of p.d.f. (with different values of
$Sk$) with the four types of PS (two values of each of
$p$ and
$\lambda _0$), 12 different roughness topographies are studied in the present work, which are summarized in table 1. Selected patches of all 12 roughness topographies on surfaces with size
$2.4H \times 0.8H$ are displayed in figure 3. Above each roughness map, the one-dimensional (1-D) roughness profile at
$z=0.4H$ along the streamwise direction is shown.
Roughness maps with configuration
$M1-500$: (a–d) negatively skewed, (e–h) zero skewness and (i–l) positively skewed; for (a,c,e,g,i,k)
$p=-1$ and (b,d,f,h,j,l)
$p=-2$. Colour indicates height. The 1-D roughness profile at
$z=0.4H$ (red dashed line) is shown above each roughness map.

Roughness topographical properties: skewness
$Sk=(1/k_{{rms}}^3)\int _S(k-k_{{md}})^3\,\mathrm {d}S$, effective slope
${ES}=(1/S)\int _S|\partial k/\partial x|\,\mathrm {d}S$, mean absolute height
$k_{a}=(1/S)\int _S|k-k_{{md}}|\,\mathrm {d}S$, r.m.s. height
$k_{{rms}}=\sqrt {(1/S)\int _S(k-k_{{md}})^2\,\mathrm {d}S}$, and meltdown height
$k_{{md}}=(1/S)\int _Sk\,\mathrm {d}S$ measured from the deepest trough, where
$S$ is the wall-projected surface area. The values of
${\rm \Delta} U^+$ and
$d/k$ are computed in full channels at
$Re_\tau =500$.

For each roughness sample, simulations in full-span and minimal channels are carried out. For minimal channels, the spanwise size
$L_z$ is the main subject of the study. As the log-layer flow structures are set by the spanwise dimension
$L_z$ (Flores & Jiménez Reference Flores and Jiménez2010), that is often most critical in terms of reducing the computational cost. The spanwise sizes of the present minimal channels are designed to fulfil the three criteria set by inequalities (1.1a,b). The first criterion (
$L_z^+>100$) stems from the fact that the computational domain must accommodate the near-wall cycle of turbulence and, unlike the other two criteria, is independent of the roughness topography. The second criterion (
$L_z^+\geq \tilde {k}^+/0.4$) ensures that roughness can be included in the ‘healthy turbulent’ zone under the critical height. We translate the criteria for a realistic roughness by replacing
$\tilde {k}$ for the sinusoidal roughness by the characteristic roughness height
$k$ for any arbitrary roughness. The third criterion states that the minimal channel should contain the roughness wavelength. However, for realistic surfaces, a single characteristic wavelength is not naturally determined. A conservative choice for this limit of the channel width can be the largest in-plane length scale, which is
$\lambda _0$ in the present study. This ensures that all wavelengths present in the roughness topography are included in the spanwise domain. Recalling the aim of reducing the cost of the roughness simulation, however, we seek a less conservative choice, in which some of the larger wavelengths are excluded. Particularly, for simulations of engineering roughness, it is often impracticable to include extremely large roughness scales. To formalize our choice, we denote the largest spanwise wavelength that a domain can accommodate as
$\lambda ^*$, and calculate the portion of surface energy that larger wavelengths contribute to the original roughness as
\begin{equation} \varPhi_c\left(\frac{2{\rm \pi}}{\lambda^*}\right)= \frac{\displaystyle\int_{2{\rm \pi}/\lambda^*}^{2{\rm \pi}/\lambda_1}E_k(q)\,\mathrm{d}q}{ \displaystyle\int_{2{\rm \pi}/\lambda_0}^{2{\rm \pi}/\lambda_1}E_k(q)\,\mathrm{d}q}, \end{equation}
where
$\lambda$ is the discrete wavelength. If the spanwise domain size is
$\lambda ^*$, the simulation resolves a roughness with
$\varPhi _c$ portion of the original surface variance.
In the current research, we examine a choice of spanwise channel size corresponding to half the size of the largest length scale, i.e.
$\lambda _0/2$. With the adopted power-law PS, this leads to values of
$\varPhi _c$ equal to or larger than 90 % for all samples under investigation. Hence the third criterion is replaced by
$L_z\geq \lambda _{1/2}$, where
$\lambda _{1/2}=\lambda _0/2$, which means that the simulations resolve at least 90 % of the original surface variance. The new criterion leads to a minimum channel width of
$L_z=0.8H$ for half of the investigated roughness topographies (those with
$\lambda _0=1.6H$) and
$L_z=0.4H$ for the other half (those with
$\lambda _0=0.8H$). We label the minimal channels with the former (larger) and latter (smaller) spanwise sizes as
$M$1 and
$M$2, respectively. For roughness samples with
$\lambda _0=0.8H$, simulations at both
$M$1 and
$M$2 channels are carried out. To complete the investigation, some simulations with a further reduced channel size
$M$3 (
$L_z=0.3H$) are carried out. It is worth mentioning that, since
$k=0.1H$ holds for all topographies,
$M$1,
$M2$ and
$M$3 channels satisfy the
$L_z\geq k/0.4$ criterion. For all simulation configurations, the streamwise channel size
$L_x$ is set according to inequalities (1.1a,b) once
$L_z$ is known.
Simulations are carried out at four different friction Reynolds numbers,
$Re_\tau =250$,
$500$,
$750$ and
$1000$, at fixed
$k/H=0.1$, leading to
$k^+=25$,
$50$,
$75$ and
$100$. However, the parametric study on roughness topography is conducted only at
$k^+=50$. Apart from minimal-channel simulations, conventional full-span channel simulations with the size
$L_x\times L_y\times L_z=8H\times 2H\times 4H$, labelled as
$F$, are also carried out for all roughness topographies. For the two highest Reynolds numbers, however, such large simulations are costly. Consequently, for these Reynolds numbers, the largest investigated channels are smaller than the
$F$ channel (but still larger than
$M$1,
$M$2 and
$M$3). These channels are labelled as
$M$0. Table 2 summarizes the details of all simulations carried out for rough channels. In order to provide a reference for determining the roughness function
${\rm \Delta} U^+$, additional smooth-wall simulations are also performed in
$M2$-,
$M1$- and
$F$-sized channels at
$Re_\tau =500$ (not shown in table 2). Overall, each rough-wall simulation case is defined by a combination of roughness topography and simulation configuration (channel size and Reynolds number). Throughout the article, the following naming convention is used to describe the cases:
\begin{equation} \overbrace{\underbrace{\boxed{\;G\;}}_{\text{p.d.f.}}\;\;\underbrace{\boxed{ \;2\;}}_{{-}p} \underbrace{\boxed{\;4\;}}_{\;10\lambda_{1/2}/H}}^{\text{topography}}\;\Bigg| \overbrace{\underbrace{\boxed{\;F\;}}_{\text{ch. size}}- \underbrace{{500}}_{10k^+}}^{\text{simulation configuration}} \end{equation}(i) The first character indicates the type of p.d.f.:
$G$ for Gaussian distribution,
$P$ for positively skewed (
$Sk\approx 0.48$) and
$N$ for negatively skewed (
$Sk\approx -0.48$).(ii) The second digit indicates the PS spectral slope:
$1$ for
$p=-1$ and
$2$ for
$p=-2$.(iii) The third digit represents half of the large cutoff wavelength
$\lambda _{1/2}$, with which the channel width is determined:
$4$ for
$\lambda _{1/2}=0.4H$ and
$8$ for
$\lambda _{1/2}=0.8H$.(iv) The following character(s) indicate(s) the channel spanwise size:
$F$ is for the full channel (
$L_z=4H$);
$M$1 and
$M$2 are for the larger (
$L_z=0.8H$) and smaller (
$L_z=0.4H$) minimal channels, respectively;
$M$3 utilizes the spanwise width
$L_z=0.3H$; and
$M$0 is introduced for the cases in which
$k^+=75$ and
$100$ with
$L_z=2H$ and
$1H$, respectively.(v) The last number denotes
$10k^+$, or, equivalently,
$Re_\tau$.
Summary of all simulation cases including roughness topography and simulation configurations. For all cases,
$L_y/H=2$ and
$N_y=401$. Moreover,
${*}{*}4$ indicates all roughness topographies with
$\lambda _{1/2}=0.4H$; and
$\varDelta _{y,k}^+$ indicates the grid size at the roughness height, i.e.
$y=0.1H$. Flow-through times (
$FTT=TU_b/L_x$) for the duration of statistics collection are shown in the last column, where
$T$ is the total integral time and
$U_b=(1/H_{{eff}})\int _0^{H}U\,\mathrm {d}y$ is the bulk velocity, with
$H_{{eff}}=H-k_{{md}}$.

2.4. Post-processing
The time-averaged flow field over a rough surface is heterogeneous in horizontal directions. In order to analyse the 1-D mean profile of the flow, we apply double averaging, as proposed by Finnigan & Shaw (Reference Finnigan and Shaw2008). The double-averaged velocity profile in the wall-normal direction
$\langle \bar {u}\rangle (y)$ is obtained by averaging the time-averaged velocity over wall-parallel directions, i.e.
where
$\bar {u}(x,y,z)$ is the time-averaged streamwise velocity,
$S$ is the wall-normal projected plan area (i.e. the area of the corresponding smooth wall) and the angular bracket
$\langle \cdot \rangle$ denotes horizontal averaging. The double-averaged velocity profile
$\langle \bar {u}\rangle (y)$ will be denoted as
$U$ for simplicity. The time-averaged velocity field is obtained over a long enough period of time. It is reported by Flores & Jiménez (Reference Flores and Jiménez2010) in their study of a smooth minimal channel that, due to the bursting events, the simulation time required to achieve converged flow statistics is longer than for a conventional full-span simulation. In order to achieve a converged mean velocity profile, the minimum flow-through time (
$FTT$) is chosen to be 300 for minimal channels. The mean velocity profiles are proven converged for all cases with 300
$FTT$s. Initial transients are removed from the statistical integration.
The influence of roughness on the mean flow can be accounted for by a modified coefficient of viscosity
$\nu _e$ beyond the region where the shape of the velocity profile is affected by roughness, i.e. outer layer (Perry & Joubert Reference Perry and Joubert1963). This
$\nu _e$ can be interpreted as a downward shift in the logarithmic layer in the inner-scaled streamwise velocity profile relative to the smooth case. This downward velocity shift in the logarithmic region is referred to as roughness function
${\rm \Delta} U^+$ (Hama Reference Hama1954; Clauser Reference Clauser1956), which is further confirmed by a number of roughness studies (e.g. Schultz & Flack Reference Schultz and Flack2009). In other words, as a result of the outer-layer similarity (Townsend Reference Townsend1976), which states that outer-layer flow is unaffected by the near-wall events except for the effect due to the wall shear stress, the downward shift of the velocity profile is approximately a constant value in the logarithmic region and possibly beyond if the outer-flow geometry and Reynolds number are matched. Introducing the roughness function to the logarithmic law of the wall (log-law hereafter), the following can be written:
Here
$\kappa \approx 0.4$ is the von Kármán constant,
$B\approx 5.2$ is the log-law intercept for the smooth wall,
$d$ indicates the zero-plane displacement, which will be discussed in detail in the following section, and the superscript
$+$ indicates scaling in wall units. Based on the pioneering work by Nikuradse (Reference Nikuradse1933), the roughness function in the fully rough regime is solely a function of the inner-scaled roughness size
$k_s^+$ for the sand-grain roughness according to
Equation (2.9) is the basis for the definition of ‘equivalent’ sand-grain roughness (also denoted by
$k_s$) for an arbitrary roughness with the same roughness function.
In the present work, the roughness function
${\rm \Delta} U^+$ is calculated as the mean offset of the inner-scaled mean velocity profile over the logarithmic layer for the cases with
$Re_\tau \approx 500$. Since corresponding smooth channels in
$M2$,
$M1$ and
$F$ with matched
$Re_\tau \approx 500$ are available, this quantity is calculated by direct comparison to the smooth case. For those cases with varied
$Re_\tau$, further smooth-channel simulations with matched
$Re_\tau$ are required in order to derive the corresponding profile shift, which causes unfavourable computational effort. Having in mind that the log-law applies for minimal smooth channels under critical height
$y_c$, a good approximation of the velocity profile in the log region of smooth channels can be drawn from the log-law; thus
${\rm \Delta} U^+$ is estimated by the velocity shift at the critical height
$y_c^+=0.4 L_z^+$ relative to the log-law
$U^+=(1/\kappa )\ln (y^+)+5.2$, where
$\kappa =0.4$.
Finally, it must be noted that, unlike a smooth channel, the origin of the wall-normal coordinate for the log-law is not naturally defined for a rough wall. In this regard, Jackson (Reference Jackson1981) suggested the use of the moment centroid of the drag profile on a rough surface as the virtual origin for the logarithmic velocity profile. The definition of the virtual wall zero-plane displacement
$d$ in the present work follows Jackson's method.
3. Results
3.1. Evaluation of minimal-channel simulations
As reviewed in the previous section, while a few of the previous studies applied the minimal-channel approach for irregular roughness, they did not systematically examine it for this type of roughness. For instance, in the work by Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), a validation case is carried out by duplicating the irregular roughness structure in both wall-parallel directions. With such an approach, the minimal-channel effect considers only the effect of fluid-domain expansion or reduction but might lack the roughness-related effect due to the repetitiveness of the rough surface. In order to comprehensively assess the applicability and limits of the minimal channel, first the simulations of the roughness topography
$G24$ with variation of channel size at matched
$k^+\approx 50$ are discussed in § 3.1.1. This is followed by the results for all different roughness topographies in § 3.1.2. As mentioned before, the roughness generation process in the present research has a random nature, where only statistical properties are prescribed. To understand if a mere ‘statistical representation’ can lead to a unique flow response, simulations are carried out for eight random realizations of the roughness topography
$G24$ with
$k^+\approx 50$. In studying different random realizations with identical p.d.f. and PS, we believe we could have contributed to answering this fundamental question. The results are presented in § 3.1.3. Finally, in § 3.1.4, one roughness topography is studied in a wide range of roughness Reynolds numbers (
$k^+=25\text {--}100$) in order to assess the prediction of the minimal channel in different rough regimes. Here the roughness topography
$G24$ is evaluated using minimal channels
$M$2 and
$M$1 and (pseudo-)full-span channels
$F$ and
$M$0.
3.1.1. Minimal channels with different spanwise sizes
The 3-D schematic representations of minimal channels
$M2$ and
$M1$ as well as the full-span channel
$F$ used for simulations at
$k^+\approx 50$ are shown in figure 4(a);
$M3$ is not shown for simplicity. The hatched pattern indicates where the roughness is mounted. Roughness topography
$G24$ in the full-size simulation is shown in figure 4(b). For a direct comparison, the boundaries of minimal channels
$M1$ and
$M2$ are represented by the black frames. The pseudo-random surfaces for each configuration are generated independently. That is, for a specific topography, the surface height map in each simulation is unique, but they all share identical statistical properties.
Comparison of channel sizes. (a) Schematic simulation domain of
$F-500$,
$M1-500$ and
$M2-500$; the hatch pattern represents roughness. (b) Roughness map of
$G24F-500$; black frames indicate minimal channels
$M1-500$ and
$M2-500$.

The inner-scaled velocity profiles obtained from roughness topography configuration
$G24$ with
$k^+\approx 50$ are shown in figure 5(a). The coloured dashed lines are the velocity profiles extracted from the smooth channel. The colour indicates the size of the channel. The critical heights of minimal channels
$M2-500$ and
$M1-500$, i.e.
$y_c^+=0.4 L_z^+=80$ and 160, are illustrated by black vertical dashed lines in figure 5(a); while the critical height of
$M3-500$, i.e.
$y_c^+=60$, is shown with a grey vertical dashed line. It can be observed from the figure that the minimal-channel cases
$M2-500$ and
$M1-500$ successfully reproduce the velocity profile of a conventional full-span channel under the critical height
$y_c$. The velocity profiles deviate above the critical height
$y_c$ due to the nature of minimal channels. For
$G24M3-500$, however, some discrepancy of the profile can be observed even under its critical height.
Simulation results of roughness type
$G24$: (a) mean velocity profiles (——, rough; - - - -, smooth), and (b) velocity offset profiles. The red vertical dash-dotted line indicates roughness height measures from the zero-plane displacement
$(k-d)^+$.

In figure 5(b), the velocity offset profiles for
$F$,
$M1$ and
$M2$ channels are displayed. The velocity offset profiles are obtained by subtracting the rough-channel velocity profile from each corresponding smooth-channel velocity profiles. The velocity offset profile from minimal channels
$M1$ and
$M2$ and full-span channel in figure 5(b) show excellent agreement. Consequently, it seems like the velocity offset is not meaningfully affected by the absence of large wavelengths in the spanwise direction with a small contribution to the roughness height PS. However, this does not hold for channel
$M3$, where
$\varPhi _c(2{\rm \pi} /L_z)=82\,\%$. A similar parameter study performed for roughness topography
$G28$ (not shown here) revealed that the velocity offset starts to deviate for channel
$M2$ (
$\varPhi _c(2{\rm \pi} /L_z)=84\,\%$). In both cases, the deviation of the velocity profiles starts when the contribution of excluded large wavelengths in the roughness height spectrum is larger than 10 %. In contrast to the minimal domain guidelines (1.1a,b) that were based on the single sinusoidal wavelength, the generalization of this idea is not entirely clear for irregular roughness, which contains a wide range of wavelength. Following the findings that the large, undulating scales do not contribute to drag (Barros et al. Reference Barros, Schultz and Flack2018), the present paper demonstrates that an a priori rule of thumb is 90 % of the surface variance. This is checked a posteriori by comparing the coherence spectra between surface variance and drag to show that the drag-carrying physics are resolved, which will be discussed in § 3.2.
In the present work,
${\rm \Delta} U^+$ of minimal channels are obtained by averaging the mean velocity offset from each critical height
$y_c^+$ to half of the channel half-height
$0.5H^+$, while the full-span channel is averaged in the region
$y+=80\text {--}250$. This gives
${\rm \Delta} U^+=6.3$ for case
$G24F-500$ and
${\rm \Delta} U^+=6.4$ and
$6.2$ for cases
$G24M1-500$ and
$G24M2-500$, respectively.
3.1.2. Minimal channels for different roughness topographies
Applying the same analysis to all topographies at
$k^+\approx 50$, the roughness functions
${\rm \Delta} U^+$ are calculated and represented in figure 6(a). In this figure,
${\rm \Delta} U^+$ predicted by the minimal channels are compared with the prediction by the full-span channels with matched topographical property. It should be mentioned that, for the case
$G24M2-500$, multiple simulations are carried out for the purpose that will be discussed in § 3.1.3. Therefore, the roughness function of
$G24M2-500$ is the mean roughness function
$\overline {{\rm \Delta} U^+}$ over the
$G24M2-500$ simulations. In figure 6(a), the
${\pm }5\,\%$ disagreement interval is illustrated by the green shading around
${\rm \Delta} U_{{Mini}}^+/{\rm \Delta} U_{{Full}}^+=1$ (red line). Another key quantity widely discussed in the framework of roughness studies is the zero-plane displacement
$d$. Similar to
${\rm \Delta} U^+$,
$d$ is often used as input to roughness models, and therefore its prediction is of practical value. The predicted zero-plane displacements
$d$ in minimal channels are compared with those in full-span channels in figure 6(b). It can be observed that the minimal-channel predictions show excellent agreement with those for the conventional full-span channel, the discrepancy being less than 5 %.
(a) Roughness function and (b) zero-plane displacement predicted by minimal channels normalized with the full-span channel prediction:
$\circ$,
$M1-500$;
$\square$,
$M2-500$. The green shading indicates the prediction error interval of
${\pm }5\,\%$. In panel (a), the
$99\,\%$ confidence interval for the different simulations of case
$G24M2-500$ is shown as an error bar.

Consistent predictions of the roughness function
${\rm \Delta} U^+$ indicate the capability of the minimal channels to reproduce the roughness function
${\rm \Delta} U^+$ of the irregular pseudo-realistic roughness even if a certain range of larger wavelengths are excluded. Obviously, in minimal-channel simulations, large turbulent structures in the outer layer cannot be resolved; therefore, a non-physical wake behaviour is observed in the outer layer of the minimal-channel mean velocity profile (see figure 5). However, as the present results suggest, capturing the near-wall turbulence in the minimal channel is adequate for the prediction of roughness function
${\rm \Delta} U^+$ as a quantification of skin-friction drag. This, obviously, cannot be generalized to all aspects of the turbulent flow.
3.1.3. Effect of randomness
In the present work, roughness is generated following a pseudo-random process with prescribed p.d.f. and PS. As a result, individually generated rough surfaces with identical statistics are not deterministically identical. This randomness can be a source of uncertainty when pseudo-random roughness is used as a surrogate of realistic roughness (and possibly explaining the scatter observed in figure 6). The pseudo-random roughness generation process can also be considered an imitation of random roughness formation processes in nature or industry. Hence it can be used to shed light on whether a statistical representation of stochastic roughness is adequate to predict the flow response.
To this end, eight rough surfaces corresponding to
$G24$ topography are generated independently. The realization of each randomly generated surface is unique, while the statistical properties are virtually identical. The averaged value of the roughness function over the eight samples calculated at
$k^+\approx 50$ in minimal channel
$M$2 is
$\overline {{\rm \Delta} U^+}_{G24M2-500}=6.43$, while the 99 % uncertainty interval of all values is 0.31. This averaged value is shown in figure 6, along with the uncertainty interval. One can observe that the uncertainty bar well encompasses the
${\rm \Delta} U^+$ in the full channel. This can be taken as an indication that the minimal-channel prediction can approximately converge to the exact value if the main uncertainty due to randomness is ruled out. Nevertheless, as stated before, the error associated with one random realization is still considerably low. Additionally, one cannot rule out a minor influence due to other factors, e.g. the nature of turbulence in the minimal channel, but the present data suggest those influences to be minor, if present. It is observed in figure 6 that the 99 % confidence bar lies in the green shaded area – the 5 % error range. Overall, it can be stated that DNS in minimal channels with matched roughness statistics is an accurate tool for the prediction of
${\rm \Delta} U^+$ of realistic roughness, apart from a small discrepancy, which is arguably linked to the effect of randomness.
It is appropriate at this point to recall that the p.d.f. and PS are effectively ‘reduced-order’ representations of the actual roughness geometry. The current results suggest that this representation can almost uniquely determine the dynamic response of flow for the studied type of roughness topographies as far as
${\rm \Delta} U^+$ is concerned. However, since a reduced model does not contain all information, minor differences in
${\rm \Delta} U^+$ among surfaces with the same p.d.f. and PS is not an unexpected observation. In other words, while the generated roughness is controlled in a global sense by the p.d.f. and PS, the local distribution of its features can be affected by randomness. For example, the occurrence of clustered or streamwise-aligned roughness peaks can lead to attenuation of drag due to the sheltering effect, which will be discussed in § 3.2.
3.1.4. Minimal channels in transitionally and fully rough regimes
The values of
${\rm \Delta} U^+$ reported for the simulations with
$k^+\approx 50$ suggest that the flow probably lies in the border between the transitionally and fully rough regimes. In order to ensure that minimal channels deliver acceptable predictions in a wide range of scenarios, including both regimes, in this section we study one roughness topography (
$G$24) in a range of roughness Reynolds numbers
$k^+\approx 25\text {--}100$. Both minimal-channel simulations
$M2$ and
$M1$ as well as large-span channel simulations
$F$ and
$M0$ are carried out. The simulation set-ups are summarized in table 2. Mean velocity profiles are shown in figure 7(a), while roughness functions
${\rm \Delta} U^+$ plotted against
$k_s^+$ are shown in figure 7(b). The inner-scaled equivalent roughness height
$k_s^+$ on the abscissa of the latter is obtained by scaling the calculated
$k_s$ (as explained below) with the viscous length scale
$\delta _\nu$ at different
$Re_\tau$. One can observe from figure 7(a) that each velocity profile deviates above the respective critical height
$y_c^+=0.4 L_z^+$, which is not shown for clarity. Based on these velocity profiles, roughness functions are obtained by calculating the velocity difference at each critical height relative to the log-law
$U^+=(1/0.4)\log (y^+)+5.2$. The equivalent sand-grain roughness
$k_s$ is calculated by fitting the roughness function to the asymptotic roughness function in the fully rough regime of Nikuradse sand-grain roughness as shown in figure 7(b). In doing so, we obtain an equivalent sand-grain roughness size of
$k_s\approx 1.05k$ for roughness topography
$G24$. Notably, the calculated values of the roughness function from both minimal and full channels show an excellent agreement. Furthermore, it can be observed that
${\rm \Delta} U^+$ asymptotically approaches the fully rough regime at
${\rm \Delta} U^+\approx 6$ for both minimal and full channels.
Simulations at different
$k^+$ from 25 to 100 for fixed
$k/H$. (a) Mean velocity profiles; line colour gradually changes from grey to black with increasing
$Re_\tau$; ——,
$F$ and
$M0$; - - - -,
$M1$;
$\cdot \cdot \cdot \cdot \cdot \cdot$,
$M2$. (b) Roughness function. Data from Nikuradse's (Reference Nikuradse1933) uniform sand-grain roughness and the Colebrook relation provided for industrial pipes (Moody Reference Moody1944) are added for comparison.

3.2. Roughness surface force
In the present simulations, the IBM introduces a volume force within the solid area imposing zero velocity and hence represents the action of pressure and viscous drag force. One of the advantages of the IBM is the explicit representation of the localized hydrodynamic force exerted by roughness (Chan-Braun, García-Villalba & Uhlmann Reference Chan-Braun, García-Villalba and Uhlmann2011), here referred to as ‘surface force’. In this section, we investigate the link between the mean surface force distribution and the roughness height distribution. Given the satisfactory performance of the minimal channel demonstrated in § 3.1, the following analysis is carried out based on the results achieved from minimal channels. Previous studies on irregular roughness report that a certain range of roughness scales is dominant in generation of skin friction (Barros et al. Reference Barros, Schultz and Flack2018). A deeper insight into the contribution of different roughness scales to the drag force is the aim of this section.
To this end, the local forcing map
$\boldsymbol {f}(x,z)$ is obtained by time-averaging the IBM force field
$\boldsymbol {f}_{{IBM}}(x,y,z,t)$ and integrating the force in the wall-normal direction
$y$:
Here
$\boldsymbol {f}(x,z)=(f_x,f_y,f_z)^{\rm T}$ is the force vector and
$f_x$,
$f_y$ and
$f_z$ are the streamwise, wall-normal and spanwise force components, respectively. One should note that, precisely speaking,
$\boldsymbol {f}$ equals force per unit density and wall-projected area. Nevertheless, as we are interested in its trend rather than its absolute value, this quantity always appears in a normalized form; hence, for brevity we refer to it as ‘force’. The visualization of the normalized forcing map for all
$M2-500$ cases with their roughness distribution maps are shown in figure 8. The entire set of surface force distributions show spanwise-elongated coherent areas of negative forcing.
(a) Roughness and (b) surface force distribution pairs in
$M$2. The exposed surface derived from the 1-D sheltering model is marked by red contour lines on the roughness distribution maps in panel (a). The separation lengths of the force peaks obtained from the autocorrelation are represented by the red bars in the upper left corner of the surface force maps in panel (b).

Furthermore, force distributions in the streamwise direction at
$z=0.2H$ are displayed for three cases,
$N24M2-500$ (negative skewness),
$G24M2-500$ (Gaussian) and
$P24M2-500$ (positive skewness), in figure 9 along with the surface height functions at the same location. In this figure, the solid blue line represents the normalized negative force profile,
${-f_x}/{f_{x,{rms}}}$, while the dashed red line represents the corresponding normalized roughness profile,
${k}/{k_{{rms}}}$. Comparing different roughness topographies, it can be observed that the Gaussian surface demonstrates a larger number of extreme force peaks than the surfaces with asymmetric p.d.f. As expected, the peaks in the surface force are mostly collocated with the peaks in roughness height. A sudden rise in the negative force is expected when the mean flow impinges on the windward side of the roughness element followed by a rapid drop on the leeward side, which is also observed in figure 9. Interestingly, the force peaks are much narrower than the surface height peaks, which can be attributed to the separation of flow behind the roughness peak. Another notable observation is that the pronounced force peaks show a much longer streamwise separation than the height peaks. Such an observation can be linked to the sheltering effect, which causes a significant reduction in the flow momentum in the wake of a tall roughness element.
Normalized force
$f_x$ and roughness distribution profile for
$N$24
$M$2 (a),
$G$24
$M$2 (b) and
$P$24
$M$2 (c) at
$z=0.2H$. The sheltering effect modelled by Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016) is illustrated by grey dashed lines, and the unsheltered surfaces are marked by the green profile.

Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016) investigated the sheltering effect on a surface roughened by rectangular-prism roughness elements and argued that, once the region sheltered by the upstream roughness element covers the neighbouring elements, the surface drag decreases. They suggested that an attenuation parameter for skin friction should incorporate the ‘shadowed area’. To provide further insight into the present observations, we apply the wake expansion model proposed by these authors – with some simplification – to the roughness profiles in figure 9. According to Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016), the streamwise slope of the sheltered region from the separation point down to the ground can be calculated from the wake expansion rate by
$\tan \theta =C_\theta u_\tau /U_h$, where
$u_\tau$ is the friction velocity,
$U_h$ is the velocity at roughness element height,
$C_\theta =1-\tfrac {2}{3}(1-h/w)$ is the shape parameter of the roughness and
$h/w$ denotes the aspect ratio of the rectangular prisms. With the aim to investigate the underlying physical mechanism of the sheltering model, we expand the use of the model – which is obtained based on rectangular-prism roughness – to more realistic irregular roughness. To approximate the expansion rate of irregular roughness in the present work, we use the double-averaged mean velocity at each roughness peak height as
$U_h$. We also replace
$h$ with the characteristic roughness height
$k_{99}=0.1H$, and replace
$w$ with the spanwise integral length scale of the roughness
$L_{k,z}\approx 0.05H$, which will be defined in the following section. Using these values,
$C_\theta =1.7$ is obtained and the resulting shadowed area in figure 9 is indicated by the grey dashed lines. In the same figure, the ‘exposed’ areas on the roughness peaks are highlighted by green lines. It is clear that the extreme force peaks coincide almost exclusively with the green areas, and the shadowed areas rarely produce any significant local force. This can be an indication of the applicability of the wake expansion model to irregular roughness. The 1-D sheltering model is applied to the 2-D roughness distribution in figure 8. The exposed roughness surface is outlined by the red contour lines. One can observe that the exposed surface contours match well with the localization of the surface force. Notably, the spanwise-elongated patterns of the surface force and their streamwise separation can be well reproduced with the help of the sheltering model. This finding can also be an indication of the feasibility to predict the local drag with a knowledge of roughness structure, which has a predictive potential for more complex problems, e.g. inhomogeneous and anisotropic roughness structures.
To shed further light on the surface force patterns and sheltering effect, the streamwise autocorrelation functions of the streamwise surface force
$\bar {R}_{ff}(\varDelta _x)$ for the
$N24M2$,
$G24M2$ and
$P24M2$ cases are shown in figure 10. The streamwise autocorrelation function of the surface force is defined as
A rapid drop of autocorrelation in the vicinity of zero separation – a result of the narrow peaks in the surface force distribution – is observed in figure 10. Additionally, a mild but clear positive peak in the autocorrelation function, as marked by filled red circles, is observed at a separation of approximately (0.3–0.6)
$H$. This value is likely to be related to the streamwise distance between the force peaks. The locations of the second autocorrelation peaks for the rest of the studied cases are marked by open red circles in figure 10 without showing the autocorrelation functions for better clarity. For visual comparison, we also indicate these values by red bars on the upper left corner of the respective surface force maps in figure 8(b). Here it can be confirmed that, roughly speaking, the lengths of the bars are similar to the separation between the spanwise-elongated areas with high surface force.
Autocorrelation function of the streamwise surface force component as a function of streamwise separation
$\varDelta _x$. The first correlation peaks of the three autocorrelation functions are marked by red filled points, whereas the peaks of the remaining cases studied are marked by open red circles.

3.2.1. Correlation between surface force and roughness height
In this section, the link between the streamwise component of the surface force,
$f_x(x,z)$, and the roughness height distribution,
$k(x,z)$, is analysed by means of the correlation function between the two quantities. The correlation function
$R_{kf}(\varDelta _{x})$ is calculated along the streamwise direction, followed by averaging in the spanwise direction:
where the subscript
$rms$ and the overbar indicate the r.m.s. value and the mean value, respectively.
The calculated values of the correlation coefficient
$\bar {R}_{kf}(\varDelta _x=0)$ for the studied topographies are summarized in table 3. A negative sign of the correlation indicates that high roughness peaks are correlated with negative surface force (force directed against the streamwise mean flow), as expected. It is observed that the negatively skewed topographies show lower correlation coefficients, which can be linked to the fact that this type of roughness is rather prone to generation of recirculation and separation zones in the surface valleys and indentations, so the responding force is less localized in those areas. In contrast, for the positively skewed topographies, higher correlation coefficients are observed due to the peak-dominated structures, in which the protruding parts of roughness are directly responsible for the generation of a localized drag force.
Cross-correlation coefficients
$\bar {R}_{kf}(\varDelta =0)$ and
$\bar {R}_{kf,{exp}}(\varDelta =0)$.

The correlation coefficients of the exposed surface with the surface force distribution,
$\bar {R}_{kf,{exp}}$, are also estimated, and are shown in table 3. Hereby only the exposed surfaces, i.e. the surface areas that are marked by red contour lines in figure 8, are kept on the height map, while the sheltered surfaces are replaced by zero elevation. A noticeable increase in the correlation coefficient can be observed for all cases, especially for negatively skewed roughness, where the correlation is increased by approximately 35 %, whereas the increases for Gaussian and positively skewed roughness are 20 % and 16 %, respectively. This, once again, highlights the importance of the sheltering effect and exposed roughness areas for the generation of surface force.
Furthermore, we calculate the integral length scales of the streamwise surface force (
$L_f$) and the roughness height distribution (
$L_k$). The integral length scales are calculated in a similar way to that proposed by Quadrio & Luchini (Reference Quadrio and Luchini2003) using the following expression for the integral length scale of roughness height:
Here
$\varDelta$ is the separation in either the streamwise or spanwise direction and
$L^{{corr}}_{0.2}$ is the separation at which the autocorrelation function drops below the arbitrary value of
$0.2$. The integral length scale for surface force
$L_f$ is computed in a similar fashion.
The calculated values can be regarded as a scale for the width of the roughness elements or force peaks. Both integral length scales are calculated for different cases and plotted in figure 11 as functions of topographical properties, i.e. PS slope
$p$ and
$\lambda _0$. The streamwise integral length scales are plotted in figure 11(a–c), while the spanwise integral length scales are plotted in figure 11(d–f), grouped by
$Sk$. Square symbols represent
$L_k$ while circles represent
$L_f$. It is observed that the surface force has a smaller streamwise integral length scale than the surface height, which is in line with the qualitative observation of the very narrow force peaks in figure 9. As expected, topographical parameters show a clear impact on the integral length scales of roughness height. In contrast to this, the streamwise integral length scale of force
$L_{f,x}$ does not show a strong sensitivity to the considered roughness variation. A notable observation is that the spanwise integral length scale of force
$L_{f,z}$ is more sensitive to the topographical changes than the streamwise length scale. The value of
$L_{f,z}$ is comparable to the surface height integral length scale
$L_{k,z}$. Based on the limited data points in our dataset, the two quantities show similar trends with
$p$ and (to some extent)
$\lambda _0$: a higher value of
$L_{f,z}$ is obtained for the roughness with
$p=-2$ and
$\lambda _0=1.6H$. Unlike the isotropic behaviour of the roughness height function, clearly illustrated by the comparable streamwise and spanwise integral length scales, the distribution of surface force is observed to be strongly anisotropic. The fact that the integral length scale of surface force is different in the
$x$- and
$z$-directions is the quantitative manifestation of spanwise-elongated coherent areas of surface force observed in figure 8.
Integral length scales
$L_k$ and
$L_f$ as functions of
$p$, grouped by
$Sk$: (a–c) streamwise integral length, and (d–f) spanwise integral length. Left axis, normalized by
$H$; right axis, normalized by
$k$. The squares represent
$L_k$ while circles represent
$L_f$. Black,
$\lambda _0=0.8H$; grey,
$\lambda _0=1.6H$.

3.2.2. Coherence function of surface force and roughness
To further understand the correlation between the force and height distributions at different roughness length scales, the coherence function between the two distributions is calculated. The coherence function represents the correlation of the force distribution with the roughness height distribution as a function of wavenumber,
where
$E_{kf}(q)$ represents the cross-PS of roughness topography
$k(x,z)$ and force map
$f_x(x,z)$, while
$E_{f}(q)$ represents the PS of
$f_x(x,z)$. Power spectra are calculated based on 1-D distribution profiles along the streamwise and spanwise directions, and the mean coherence function
$\overline {\gamma _{kf}^2}$ is obtained by averaging each parallel signal pair.
Figure 12 shows the mean coherence function of Gaussian surfaces in (a) streamwise direction and (b) spanwise direction as a function of wavelength
$\lambda ={2{\rm \pi} }/{q}$; the upper axis shows the inner-scaled wavelength
$\lambda ^+$ at
$Re_\tau =500$. It is worth reminding that, for all topographies, the smallest in-plane roughness scale is prescribed to be
$\lambda _1=0.08H$. This is clearly related to the observation that coherence functions are significantly smaller below
$\lambda \approx 0.08H$ (
$\lambda ^+\approx 40$). Above this threshold, coherence functions increase and retain high values until a certain wavelength, which is roughly at
$\lambda \approx (0.3\text {--}0.6)H$ (
$\lambda ^+\approx 150\text {--}300$) for the studied cases. With further evolution of the coherence function to larger wavelengths, the coherence decreases monotonically. Similar observations are made for negatively and positively skewed roughness. The force becomes less correlated with the surface features at very large scale, or, in other words, very large wavelengths in the streamwise direction do not contribute to the generation of surface force. These length scales might be related to the ‘wavy roughness’ concept stemming from the observations by Schultz & Flack (Reference Schultz and Flack2009). Barros et al. (Reference Barros, Schultz and Flack2018) stated that these length scales can be filtered out in regard to determining the skin friction.
Mean coherence function
$\overline {\gamma _{kf}^2}$ as a function of
$\lambda =2{\rm \pi} /q$ normalized by
$H$: (a) streamwise and (b) spanwise. The value
$\lambda _1=0.08H$ is marked by vertical dashed lines on the left and right side, respectively; and
$\lambda ^+$ for the current cases
$Re_\tau \approx 500$ is shown on the upper axis. The length scale of force separation obtained from its autocorrelation functions are marked by red circles on each coherence function, respectively.

Notably, the streamwise wavelength at which the coherence starts to drop has a similar value to the streamwise separation distance of the surface force peaks. A comparison between the coherence dropping wavelength
$\lambda _{{Coh}}$ and the length scale of force peak separation
$\lambda _{f}$ is conducted in figure 13. It can be observed that all data points are clustering around
$\lambda _{{Coh}}=\lambda _{f}$ (dotted line), indicating a clear correspondence of these length scales. The significant coherence at relatively small length scales might be linked to the interaction of roughness structure and sheltering. As discussed before, the occurrence of extreme force peaks is strongly determined by the roughness areas that are exposed in the high-momentum flow, or outside the sheltering. Thus, streamwise recurring force peaks caused by sheltering can be found in figure 8. Less prominent force peaks can be found between two successive extreme peaks, which contributes to the coherence function at small wavelength. Furthermore, the roughness structures whose length scales are comparable to the distance between two successive extreme peaks, i.e.
$\lambda _f$, show significance in the coherence function at corresponding wavelength. Beyond this length scale, no larger force peak separation can be found, and the coherence function keeps decreasing into the long-wavelength region. Figure 12(b) demonstrates that, unlike the tortuous behaviour exhibited by the streamwise mean coherence function, the spanwise mean coherence function shows a monotonically increasing trend before a plateau at larger wavelengths.
The length scale of force peak separation detected from the autocorrelation functions
$\lambda _{f}$ compared with the coherence dropping wavelength
$\lambda _{{Coh}}$.

3.3. Effect of roughness topographical properties
In § 3.1, we discussed the applicability of the minimal-channel concept for characterization of realistic roughness. In doing so, we examined a relatively large number of roughness topographies, which provide a basis for studying the effect of roughness topography on the hydrodynamic properties of the surface, which is discussed in the present section. It has been already shown that p.d.f. and PS can be considered as a reduced representation of roughness topography (almost) uniquely reproducing the hydrodynamic response. The common practice in the literature is, however, to parameterize roughness in terms of a few statistical parameters – an even further reduced representation. Aiming to establish a link with the existing literature, we adopt and examine this approach in § 3.3. We first show the trends of
${\rm \Delta} U^+$ and zero-plane displacement with some statistical parameters in § 3.3.1. Then some existing roughness correlations are assessed based on the present results in § 3.3.2.
3.3.1. Effect of statistical parameters on
${\rm \Delta} U^+$ and
$d$
An overview of the roughness function
${\rm \Delta} U^+$ for all topographies (see table 1) is plotted in figure 14(a), where the roughness function
${\rm \Delta} U^+$ is shown as a function of two of the investigated roughness parameters, i.e. skewness
$Sk$ and PS slope
$p$. As investigated by Flack et al. (Reference Flack, Schultz and Barros2020), positively skewed rough surfaces give higher skin friction than non-skewed or negatively skewed roughness. In general,
${\rm \Delta} U^+$ reaches a higher value with
$p=-1$. As illustrated in figure 2(b), at
$p=-2$ larger wavelengths contribute more to the roughness and vice versa at
$p=-1$. These findings agree with the study by Barros et al. (Reference Barros, Schultz and Flack2018) and the results in § 3.2 highlighting that larger horizontal length scales contribute less to the hydrodynamic drag. It should be noted that the effective slope of roughness is larger for
$p=-1$ compared to the corresponding cases with
$p=-2$. Furthermore, it can be observed that the surfaces with
$\lambda _0=1.6H$ show stronger sensitivity to the change of PS slope
$p$ than those with
$\lambda _0=0.8H$.
(a) Effect of roughness topographical properties on the
${\rm \Delta} U^+$ prediction from minimal channels: black,
$\lambda _0=0.8H$; grey,
$\lambda _0=1.6H$. (b) The
$d$ predictions as a function of
$ES$, grouped by
$Sk$.

Although the observed trends with selected statistical parameters can, to some extent, justify the use of these parameters in predictive correlations, it is also observed that predicting skin friction based on a few roughness statistics is incomplete. Considering the roughness statistics and the values of roughness function in table 1, it is observed that different types of roughness with similar statistical properties, e.g. roughness with same values of
$Sk$ and
$k_{rms}$ in table 1, show meaningful variation of their
${\rm \Delta} U^+$ values. A better correlation for
${\rm \Delta} U^+$ is observed when
$ES$ is added to the prediction. However, inclusion of
$ES$ is still not necessarily expected to yield unique predictions. A simple illustrative example could be that roughness formed by staggered and aligned roughness elements with identical statistics can lead to significantly different values of
$k_s$, as shown in the study by Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017).
As stated before, for a rough surface, the logarithmic layer of the flow is shifted upwards with respect to the bottom plane. Thus, the origin of the wall-normal coordinate cannot be defined a priori. As a result, it is necessary to use a physically justified virtual origin for the logarithmic law of the wall. The virtual origin lies above the
$y=0$ plane at a distance equal to the zero-plane displacement
$d$. The value of the zero-plane displacement
$d/k$ following Jackson's method (Jackson Reference Jackson1981) is documented in table 1. To summarize the effect of roughness topography on the zero-plane displacement,
$d/k$ is plotted as a function of effective slope
$ES$ in figure 14(b), while the data points are grouped by
$Sk$. Even though
$ES$ is not explicitly prescribed in the present work, it is indirectly controlled by the two PS parameters
$p$ and
$\lambda _0$. It can be observed that the value of
$d/k$ increases with an increase in
$ES$ and a decrease in
$Sk$, while the skewness effect is more dominant.
3.3.2. Assessment of existing roughness correlations
In this section, results from previously introduced topographies at
$Re_\tau \approx 500$ are used to assess some of the existing roughness correlations. In fact, existing roughness correlations are developed based on a limited number of data points covering a certain region of the parameter space (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). In this section, we are particularly interested to shed light on the generalization of these correlations outside their original parameter space, which is key for a correlation to work across the wide range of rough surfaces encountered in different applications.
In the following, we assess three relatively recent correlations by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015), Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) and Flack et al. (Reference Flack, Schultz and Barros2020), each predicting
$k_s$ based on a few roughness statistical parameters. These correlations are applicable in the fully rough regime. Figure 15 visualizes the selected correlations, where the parameter space covered by the original fitting data of each correlation is represented by a red frame. In each panel the data points from the present work are depicted as symbols. Different symbol colours are used to make a distinction between the data points lying inside and outside the parameter space originally used for development of the correlation. Here the parameter space is expressed in terms of the two widely used parameters
$Sk$ and
$ES$.
Predictive correlations for
$k_s^+$. Squares indicate data points from the present work while the red frame represents the fitting data from the literature: (a) (3.6) by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015), (b) (3.7) by Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017), and (c) (3.8) by Flack et al. (Reference Flack, Schultz and Barros2020).

The required roughness statistics of the present roughness topographies are listed in table 1. First we examine the correlation proposed by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015),
which is developed based on 3-D sinusoidal roughness. In figure 15(a) the data points from the present work are located in the range of fitting data except for the topographies with
$Sk<0$. Obviously, since
$Sk$ is not used as a predictive parameter in this correlation, it returns the same predictions for different values of
$Sk$.
Furthermore, two correlations developed by Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017),
and by Flack et al. (Reference Flack, Schultz and Barros2020),
\begin{equation} \left\{\begin{array}{@{}ll} k_s=2.73k_{{rms}}(2+Sk)^{{-}0.45}, & Sk<0,\\ k_s=2.11k_{{rms}}, & Sk=0,\\ k_s=2.48k_{{rms}}(1+Sk)^{2.24}, & Sk>0, \end{array}\right. \end{equation}are examined. The original fitting data, from which these correlations are extracted, also only partly cover the current data as illustrated in figure 15.
In order to directly evaluate the models,
$k_s^+$ predicted by the correlations are normalized and plotted against the full-span DNS results in figure 16. The equivalent sand-grain sizes
$k_s^+$ of roughness from the simulations are obtained by fitting the roughness function to the fully-rough asymptote, i.e. (2.9). One should recall that the present data points can cover both transitionally and fully rough regimes, while the correlations in question are to be examined for the latter regime. Based on the result in § 3.1.4, and similar to the approach adopted by Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), an approximate value of
${\rm \Delta} U^+\approx 6$ is regarded as the threshold of the fully-rough regime and only the data points in the fully-rough regime are shown in figure 16. It has to be mentioned that this value is an approximate criterion to estimate the roughness regime. The exact
${\rm \Delta} U^+$ criterion for different types of roughness can only be achieved through comprehensive experiments.
Prediction of
$k_s^+$ compared to DNS. Red symbols indicate the present data points located within the fitting dataset. The 10 %, 20 % and 30 % error intervals are represented by green shading.

In figure 16, different error intervals in prediction of
$k_s^+$ are illustrated by various green shades. In general, a similar range of error can be observed among the predictions of all correlations in the figure where a limited number of data points lie outside the 30 %
$k_s$ error area. The fact that none of the correlations are able to perfectly reproduce the effect of topography on
$k_s$ – as already pointed out by other authors (Flack et al. Reference Flack, Schultz and Barros2020) – can be acknowledged in figure 16. Among all the correlations, those by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015) and Flack et al. (Reference Flack, Schultz and Barros2020) incorporate less geometrical information by taking one parameter related to the topography each, namely the effective slope of roughness in the former and the skewness in the latter. Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) combined both approaches. While the latter correlation delivers better predictions for some data points, an obviously superior accuracy cannot be established. One notable observation from all the correlations is that their prediction does not deteriorate particularly for data outside their original fitting range. This can be an indication that either of the models can be used with a similar level of reliability in a wider parameter space than it was originally designed for. Obviously, this statement is unlikely to hold for extreme cases outside the scope of this paper.
To summarize, this section highlights the need for a more general model, e.g. by taking advantage of advances in data-driven methods. For the future development of roughness models/correlations, it is suggested by the present work that roughness statistics that contain both the height distribution and horizontal scales of roughness need to be incorporated for generalizable predictions.
4. Conclusion
DNS is carried out for turbulent flow over irregular roughness in plane channels with reduced stream- and spanwise extents – referred to as minimal channels. Roughness topography is mathematically generated using the method proposed by Pérez-Ràfols & Almqvist (Reference Pérez-Ràfols and Almqvist2019), in which the p.d.f. and PS of the roughness map can be prescribed with high precision. Simulations are run for 12 different roughness topographies at
$k^+=50$ and for a selected topography at
$k^+=25\text {--}100$ (spanning both transitionally and fully rough regimes). For all cases, solutions are produced in full channels and one or more minimal channels. It is systematically demonstrated that the value of the roughness function for an irregular roughness with random nature can be predicted within
${\pm }5\,\%$ error using DNS in a minimal channel. This can be achieved as long as the minimal-channel dimensions follow a relaxed version of the criteria suggested for regular sinusoidal roughness by previous authors (Chung et al. Reference Chung, Chan, MacDonald, Hutchins and Ooi2015; MacDonald et al. Reference MacDonald, Chung, Hutchins, Chan, Ooi and Garcia-Mayoral2017). The relaxation concerns the condition that the channel should contain all horizontal scales of roughness. Current data suggest that accurate prediction can be achieved as long as the size of the channel is large enough to accommodate more than 90 % of the original roughness height spectral energy based on the area under the 2-D PS. This finding is particularly relevant in DNS-based characterization of realistic rough surfaces, which may contain very large wavelengths with limited contribution to the r.m.s. roughness height.
To shed more light on possible origins of the mild discrepancy between the minimal- and full-channel results, for one topography, multiple rough surfaces are generated. Owing to the random nature of the roughness generation process, these surfaces are deterministically different while being statistically identical. Simulations are carried out for these surfaces with
$k^+=50$ and slightly scattered values of the roughness function are obtained. Notably, the roughness function value for the full channel resides within the 99 % uncertainty interval of these scattered predictions. The results indicate that, at fixed p.d.f. and PS, randomness in roughness generation can lead to a small uncertainty, which is also likely to be the origin of the observed
${\pm }5\,\%$ discrepancy between the predictions of the minimal and full channels. This can be an indication that one can consider a combination of p.d.f. and PS as a reduced-order representation of roughness topography, leading to nearly unique dynamic flow responses.
In addition to global flow properties, local surface forces for different types of roughness are calculated and their correlations with respective roughness height functions are studied. It is observed that not all roughness height peaks generate force peaks. Applying the sheltering model proposed by Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016) with some assumptions, we are able to show that only ‘exposed’ (unsheltered) roughness peaks generate prominent peaks in the surface force. Notably, the spanwise-elongated patterns of the surface force and their streamwise separation can be well reproduced with the help of the sheltering model. This can be taken as a clear indication of the relevance of the sheltering effect in flow over irregular roughness – e.g. for complex terrains.
To shed light on the contributions of different roughness scales to global drag, we also studied the coherence function of roughness height and surface force power spectra as a function of streamwise and spectral scales. In the streamwise direction, it was observed that coherence starts dropping beyond a certain length. This observation can be interpreted as a smaller contribution of very large roughness wavelengths to the drag force. These large roughness length scales might be related to the ‘wavy roughness’ concept stemming from the previous studies (Schultz & Flack Reference Schultz and Flack2009).
Our analysis of surface force reveals certain previously unattended facts about roughness-induced drag, e.g. reduced coherence between friction and roughness height at large scales. Notably, the wavelength at which the coherence starts dropping is shown to be related to the separation between the peaks of surface force, which is linked to the sheltering effect itself. Unlike the streamwise direction, the coherence function does not drop in the spanwise direction for the cases studied in this paper.
As stated above, the present results suggest that an accurate yet computationally economical framework for the characterization of irregular, realistic rough surfaces is in hand. Such a framework can, for example, be a basis for generation of the large databases required for future ‘data-driven’ roughness correlations. While this can be considered an obvious future research direction, in the present paper, we used the results from the 12 simulated roughness topographies to study the dependence of the roughness function as well as the zero-plane displacement on some key roughness parameters. Notably, it was shown that the normalized zero-plane displacement
$d/k$ is most sensitive to the skewness of the roughness distribution (larger at smaller values of skewness), and it also mildly increases with effective slope.
Finally, we assessed a number of widely cited roughness correlations in the literature, here presented in (3.6)–(3.8). While some correlations show a certain level of success in reproducing the roughness function or equivalent sand-grain roughness compared to the DNS results (see figure 16), there is an obvious need for improvement. An interesting observation is that none of the assessed correlations show a dramatic loss of accuracy when used outside the parameter space of its original fitting data. However, even the most successful correlations can reproduce the DNS data only within
${\pm }30\,\%$ accuracy. This can arguably be the ground for a paradigm shift in the development of future roughness correlations. As mentioned before, a data-driven approach, which can account for the stochastic nature of roughness and its interaction with near-wall turbulence, may be a solution to this problem. Recently, this idea has received some attention (Brereton, Jouybari & Yuan Reference Brereton, Jouybari and Yuan2021; Jouybari et al. Reference Jouybari, Yuan, Brereton and Murillo2021) and more work in this direction is called for.
Acknowledgements
The data presented in the manuscript are openly available in the KITopen repository at doi:10.5445/IR/1000142136.
Funding
J.Y. and P.F. gratefully acknowledge financial support from the Friedrich and Elisabeth Boysen Foundation (BOY-151). This work was performed on the supercomputer ForHLR Phase 2 and the storage facility LSDF funded by the Ministry of Science, Research and the Arts, Baden-Württemberg, and by the Federal Ministry of Education and Research.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Visualization of grid resolution
In the present work, the IBM is employed to introduce roughness into the simulation flow field. With the help of the IBM, the complex geometry of the roughness can be represented on a simple Cartesian grid. However, in order to completely represent the roughness with a Cartesian grid, sufficiently fine grid resolution is essential. The grid resolution in the present work is illustrated in figure 17, where one of the rough surfaces with the smallest correlation length
$L_x^{{corr}}$, i.e.
$G14M2-500$, is shown. This rough surface contains the finest roughness structure among all the cases.
Mesh visualization in the wall layer with rough surface
$G14M2-500$: (a) mesh in wall-normal direction along
$z=0.2H$; and (b) mesh in wall-parallel directions.

Appendix B. Convergence test of integration time
In order to determine the appropriate integration time for the converged mean velocity profile, mean velocity profiles over roughness
$G24M1$ at four different
$Re_\tau$ averaged with 150 and 300
$FTT$s are compared and shown in figure 18. The profiles obtained with 150
$FTT$s at four different
$Re_\tau$ collapse onto the 300
$FTT$ profiles. On the other hand, MacDonald et al. (Reference MacDonald, Chung, Hutchins, Chan, Ooi and Garcia-Mayoral2017) estimated the convergence of the mean velocity profile in minimal channels in terms of the count of captured
$y_c$-sized eddies during the simulation. The number of captured eddies can be expressed as
where
$T_{sim}$ is the total simulation time, and
$y_c=0.4L_z$ is the critical height of the minimal channels. The 95 % confidence interval of the
${\rm \Delta} U^+$ prediction is formulated as
${\rm \Delta} U^+\pm \epsilon ^+$, where
$\epsilon ^+\approx 91.4(C^*)^{-1/2}/y_c^+$. The current criteria of 300
$FTT$s corresponds to
$\epsilon ^+\approx 0.07$, 0.07, 0.05 and 0.03 for minimal channel
$M1$ at
$Re_\tau =250$, 500, 750 and 1000, respectively. Thus the integration time of minimum 300
$FTT$s in the present work is shown to be long enough for achieving a converged mean velocity profile.
Statistical convergence of minimal channel
$M1$; colour indicates different
$Re_\tau =250$, 500, 750 and 1000; lines, 300
$FTT$; squares, 150
$FTT$.

Appendix C. Definition of zero-plane displacement
The zero-plane displacement
$d$, according to Jackson (Reference Jackson1981), is placed at the centroid of the distributed drag on the roughness surface. The moment of the drag can be calculated by projecting the drag forces on a
$y$–
$z$ plane. However, it is demonstrated by Jackson that
$d$ can be calculated in terms of mean flow properties, i.e. total shear stress
$\tau _{tot}$. The time-averaged Navier–Stokes equation of the flow in the streamwise direction can be written as
Here
$\bar {u}$,
$\bar {v}$ and
$\bar {w}$ are the mean velocity components;
$\bar {p}$ is the mean pressure;
$T_{11}$,
$T_{12}$ and
$T_{13}$ are the stresses, including Reynolds stresses;
$\overline {f_{x,{IBM}}}$ is the streamwise component of the mean IBM force; and
$P_x$ is the constant pressure gradient added to the flow. With the idealized geometry proposed by Jackson (Reference Jackson1981), if (C1) is integrated over wall-parallel directions, then multiplied by
$y$ and integrated over
$y$ from the bottom
$y=0$ to the tip of the roughness
$y=k_t$, we obtain
\begin{align} & \underbrace{\int_{L_z}\int_0^{k_t}y\rho[\bar{u}\bar{u}]^{L_x}_0\,{\rm d} y\,{\rm d} z}_{{=}0}+\int_{L_z}\int_{L_x}[y\rho \bar{u}\bar{v}]^{k_t}_0\,{\rm d} x\,{\rm d} z \nonumber\\ & -\int_{L_z}\int_{L_x}\int^{k_t}_0\rho \bar{u}\bar{v}\,{\rm d} y\,{\rm d} x\,{\rm d} z +\underbrace{\int_{L_x}\int_0^{k_t}y\rho[\bar{u}\bar{w}]^{L_z}_0\,{\rm d} y\,{\rm d} x}_{{=}0} \nonumber\\ & \quad = \underbrace{\int_{L_z}\int^{k_t}_0-y[\bar{p}]^{L_x}_0\,{\rm d} y\,{\rm d} z}_{{=}0} +\underbrace{\int_{L_z} \int_0^{k_t}y[T_{11}]^{L_x}_0\,{\rm d} y\,{\rm d} z}_{{=}0} +\int_{L_z}\int_{L_x}[yT_{12}]^{k_t}_0\,{\rm d} x\,{\rm d} z \nonumber\\ &\qquad - \int_{L_z}\int_{L_x}\int^{k_t}_0T_{12}\,{\rm d} y\,{\rm d} x\,{\rm d} z+ \underbrace{\int_{L_x}\int_0^{k_t}y[T_{13}]^{L_z}_0\,{\rm d} y\,{\rm d} x}_{{=}0} \nonumber\\ &\qquad -\int_{L_z}\int_{L_x}\int^{k_t}_0yP_x\,{\rm d} y\,{\rm d} x\,{\rm d} z+ \int_{L_z}\int_{L_x}\int^{k_t}_0y\bar{f}_{x,{IBM}}\,{\rm d} y\,{\rm d} x\,{\rm d} z. \end{align}As marked in the equation, some of the terms vanish due to the periodic boundary condition in wall-parallel directions. Thus, the moment of the drag acting on the roughness can be written as follows:
\begin{align} M &={-}\int_{L_z}\int_{L_x}\int^{k_t}_0y\bar{f}_{x,{IBM}}\,{\rm d} y\,{\rm d} x\,{\rm d} z \nonumber\\ &= \int_{L_x}\int_{L_z}[yT_{12}-y\rho \bar{u}\bar{v}]_{y=k_t}\,{\rm d} z\,{\rm d} x -\int_{L_x}\int_{L_z}\int^{k_t}_0[T_{12}-\rho \bar{u}\bar{v}]\,{\rm d} y\,{\rm d} z\,{\rm d} x \nonumber\\ &\quad -\int_{L_x}\int_{L_z}\int^{k_t}_0yP_x\,{\rm d} y\,{\rm d} z\,{\rm d} x. \end{align}
Here
$M$ is the moment on the surface. Following which, the zero-plane displacement
$d=k_t-M/(\tau _wL_xL_z)$ is calculated. With the operation
$(L_xL_z)^{-1}\int _{L_x}\int _{L_z}[T_{12}-\rho \bar {u}\bar {v}]\,{\rm d} z\,{\rm d} x$, dispersive stress is included in the total shear stress and is labelled as
$\tau _{tot}$ (Jackson Reference Jackson1981). Kameda, Mochizuki & Osaka (Reference Kameda, Mochizuki and Osaka2018) calculated the displacement
$d$ by setting the wall coordinate origin at
$k_t$; thus the equation for the zero-plane displacement
$d$ can be written as
\begin{equation} d=k_t-\frac{\displaystyle\int^{k_t}_0(\tau_{tot}+yP_x)\,{{\rm d}y}}{\tau_w}. \end{equation}























































































































