2.1. Numerical method and details

We perform DNS using the open-source solver Xcompact3d (Laizet & Lamballais Reference Laizet and Lamballais2009; Laizet & Li Reference Laizet and Li2011; Bartholomew et al. Reference Bartholomew, Deskos, Frantz, Schuch, Lamballais and Laizet2020), using sixth-order compact finite differences in space combined with an explicit third-order Runge–Kutta scheme in time. We test different configurations for varying $\Lambda _s$ at two different friction Reynolds numbers ( $Re_\tau =180$ and $Re_\tau =500$ ). While the streamwise extent $L_x$ of the simulation domain is set to values that are greater or equal to those used by Neuhauser et al. (Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022) for analogous simulations, the spanwise box size $L_z$ is set alternatively to $12h$ , $8h$ or $6h$ to accomodate an integer, even number of strips for each of the tested values of $\Lambda _s$ . The $L_z=6h$ box size is preferred at high- $Re$ wherever possible to minimise the computational load of a single simulation.

Our data production pipeline consists of two stages. First, initial conditions are produced by simulating a channel flow with a spanwise roughness pattern of period $\Lambda _s$ (see figure 1 a). Rough wall sections are modelled by imposing a slip length $\ell$ for the spanwise velocity component at the wall as done by Neuhauser et al. (Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022). This results in the following Robin boundary condition:

(2.6) \begin{equation} w_{w} = \ell \, \hat {n}_{w}\cdot \left ( \nabla w \right )_{w}, \end{equation}

where the $(\cdot )_w$ subscript indicates a quantity evaluated at the wall, $\hat {n}_{w}$ a unit vector that is orthogonal to the wall and pointing into the fluid. The value of the slip length is set to $\ell ^+ = 9$ following Neuhauser et al. (Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022). After the simulation reaches a steady state, a set of $N_s$ snapshots is stored; the sample time is set to $1 h/u_p$ to ensure that snapshots are reasonably uncorrelated. As previously explained, $u_p$ is analogous to $u_\tau$ ; then, if $h$ is the maximum height of an attached eddy, its lifetime can be expected to be of the order of the eddy turnover time $1 h/u_p$ (see, e.g., Lozano-Durán & Jiménez Reference Lozano-Durán and Jiménez2014). One can then expect that the turbulent features observed in a given snapshot differ from those of the successive snapshot saved after $1 h/u_p$ .

Each of the $N_s$ saved snapshots of the secondary motions is used as the initial condition for a second simulation between smooth walls (see figure 1 b). The duration of the simulation $T_f$ is chosen to satisfactorily capture the decay of the secondary motions. Streamwise-averaged flow fields are stored every $0.01 h/u_p$ , even though a laxer time resolution could have been used in retrospect. Streamwise-averaged velocity fields are preferred to three-dimensional snapshots as the present procedure is particularly data intensive. Exploiting the several repetitions of the simulation, a total of $N_s$ fields at the same time $t$ from the initial conditions are averaged together to produce an ensemble-average of the decaying secondary motions. The whole procedure is repeated for different values of $\Lambda _s$ and $Re_\tau$ .

Table 1. Numerical details for all tested combinations of $\Lambda _s/h$ (the spanwise period of the roughness pattern) and $Re_\tau$ for our smooth (steady) and time-evolving simulations. The number of fields used to calculate statistics is indicated by $N_0=N_s$ (where $N_0$ refers to steady-state simulations, $N_s$ to time-evolving ones). $T_f$ indicates the time duration of the decaying simulation, $L_x$ and $L_z$ refer to the simulation box size in the streamwise and spanwise directions; the grid spacing is uniform in these two directions and is indicated by $\Delta x$ , $\Delta z$ respectively. The wall-normal grid spacings at the wall and centreline are instead indicated by $\Delta y_w$ and $\Delta y_c$ . The maximum in time, over each grid point and over the three spatial directions (or velocity components) of the Courant–Friedrichs–Lewy ( $\text {CFL} = V\Delta t/q$ , where $\Delta t$ is the simulation time step, $q$ is the grid spacing at some generic point in a given direction and $V$ is the velocity component at that point in the same direction) and Fourier ( $F\!o = \nu \Delta t / q^2$ ) numbers are also reported. The dot to the left of each row indicates the colour used in the following figures to indicate a given value of $\Lambda _s/h$ .

Numerical details for the complete dataset used for this study are reported in table 1. The grid spacing is normalised with the worst-case value $\delta _{v,m}$ of the viscous length scale; such a value is usually observed at the initial conditions. Additionally to the decaying simulations, two reference simulations between smooth walls ( $Re_\tau =180$ , $Re_\tau =500$ ) have been produced using the same grid used for rough simulations. In this case, $N_{0}$ indicates the number of samples used for the calculation of steady-state statistics. Notice that we adjust the number of fields ( $N_s$ , $N_0$ ) used for the computation of statistics as we change $Re_\tau$ and $\Lambda _s$ ; as for the quantification of the degree of statistical convergence, see § 3.5. Generally, $Re_\tau$ has a favourable effect on the convergence of statistics: as near-wall turbulent structures become smaller with larger $Re_\tau$ , a larger number of these features are contained in a single flow snapshot. This yields a quicker falloff of the small-scale noise. It is thus expected that a lower number of snapshots ( $N_s$ , $N_0$ ) is required for statistics to converge at high- $Re_\tau$ . The effect of the period $\Lambda _s$ depends instead on the size $L_z$ of the simulation box. As previously explained, data from the several spanwise periods contained in a single snapshot are averaged together (phase-averaging); the larger the number of periods $L_z/\Lambda _s$ in a single snapshot, the lower we expect the required $N_s$ (or $N_0$ ) to be.

3.1. Triple decomposition; momentum pathways and circulatory motions

Unlike homogeneous flows, which are best described using a Reynolds decomposition, flows featuring secondary motions are commonly described in terms of a triple decomposition. As would be done in a Reynolds decomposition, velocity fluctuations $\boldsymbol {u}^{\prime }$ are separated from the expected value $ \langle {\boldsymbol {u}} \rangle$ . Additionally, the expected value $ \langle {\boldsymbol {u}} \rangle$ is further split into its spanwise average $\boldsymbol {U} = \langle {\boldsymbol {u}} \rangle _z$ and a dispersive field $\boldsymbol {u}_d = (\tilde {u},\tilde {v},\tilde {w})$ to yield the triple decomposition. This is done as the expected value $ \langle {\boldsymbol {u}} \rangle$ of the velocity depends on the conditioned spanwise variable $\zeta$ (see § 2); using the triple decomposition, the spanwise-uniform field $\boldsymbol {U}$ is separated from the spanwise-heterogeneous dispersive field $\boldsymbol {u}_d$ ,

(3.1) \begin{equation} \left \langle {\boldsymbol {u}} \right \rangle = \boldsymbol {U} + \boldsymbol {u}_d. \end{equation}

Since only the streamwise component has a non-zero spanwise average for the present geometry,

(3.2) \begin{equation} \boldsymbol {U} = U(t,y) \,\hat {x}, \end{equation}

(3.3) \begin{equation} \tilde {u} = \left \langle {u} \right \rangle - U(t,y),\,\,\,\,\,\,\tilde {v} = \left \langle {v} \right \rangle ,\,\,\,\,\,\,\tilde {w} = \left \langle {w} \right \rangle , \end{equation}

where $\hat {x}$ is a unit vector pointing in the streamwise direction and $U$ will be referred to as the mean velocity profile. The full velocity field then reads

(3.4) \begin{equation} \boldsymbol {u}(t,x,y,z) \,=\, U(t,y)\,\hat {x} \,+\, \boldsymbol {u}_d(t,y,\zeta (z)) \,+\, \boldsymbol {u}^{\prime}(t,x,y,z). \end{equation}

Note that the averaged fields $U$ , $\boldsymbol {u}_d$ are streamwise-invariant owing to the channel geometry. This allows to further split the dispersive field $\boldsymbol {u}_d$ into two separate parts. Indeed, the continuity equation for $\boldsymbol {u}_d$ reads

(3.5) \begin{equation} \frac {\partial \tilde {v}}{\partial y} + \frac {\partial \tilde {w}}{\partial \zeta } = 0 \text {.} \end{equation}

The above equation indicates that the two-dimensional vector field given by $\tilde {v}$ and $\tilde {w}$ is divergence-less. The distribution of $\tilde {v}$ can be thus determined if $\tilde {w}$ is known (or vice versa), and the two form a single circulatory pattern. Such a cross-sectional circulatory motion is typical of secondary motions (see e.g. Neuhauser et al. Reference Neuhauser, Schäfer, Gatti and Frohnapfel2022) and will be treated separately from the remaining velocity component $\tilde {u}$ , which contains the information regarding the momentum pathways (e.g. Womack et al. Reference Womack, Volino, Meneveau and Schultz2022).

3.2. Velocity spectra of the dispersive field

In § 4, the dispersive field of the simulated flows will be inspected in real space to reveal the presence of secondary motions and their features. An additional analysis will be performed in spectral space by scrutinising the velocity spectra $\Phi _{\tilde {u}\tilde {u}}$ , $\Phi _{\tilde {v}\tilde {v}}$ , $\Phi _{\tilde {w}\tilde {w}}$ of the dispersive field. Velocity spectra are typically defined as the Fourier transform of the velocity correlation function (see, e.g., Davidson Reference Davidson2015). Such a definition cannot be used in the present case owing to the periodicity of the dispersive field: a Fourier series is used instead. As an example, the spectrum $\Phi _{\tilde {u}\tilde {u}}$ of $\tilde {u}$ is defined as

(3.6) \begin{equation} \Phi _{\tilde {u}\tilde {u}} = \frac {\mathcal {F}_z\{\tilde {u}\}^\dagger \mathcal {F}_z\{\tilde {u}\}}{\Delta \kappa _z}, \end{equation}

where $(\cdot )^\dagger$ indicates the conjugate of a complex number, $\Delta \kappa _z$ is the Fourier resolution in the spanwise direction and $\mathcal {F}_z\{\tilde {u}\}$ indicates the coefficients of the Fourier series of $\tilde {u}$ . These are defined as

(3.7) \begin{equation} \mathcal {F}_z\{\tilde {u}\} (t,y,\kappa _z) = \frac {1}{\Lambda _s} \int _0^{\Lambda _s} \tilde {u}(t,y,\zeta )\, e^{-i\kappa _z \zeta }\, \text {d}\zeta \end{equation}

where $i$ is the imaginary unit and $\kappa _z$ the spanwise wavenumber. Owing to Parseval’s theorem, the spectrum of $\tilde {u}$ can be related to the spanwise average of its energy:

(3.8) \begin{equation} \sum _{\kappa _z=-\infty }^{+\infty } \Delta \kappa _z \Phi _{\tilde {u}\tilde {u}}(t,y,\kappa _z) = \frac {1}{\Lambda _s} \int _{0}^{\Lambda _s} \tilde {u}^2(t,y,\zeta )\,\text {d}\zeta = \big \langle {\tilde {u}^2} \big \rangle _z \end{equation}

The above equation justifies the common interpretation of the spectrum as the contribution of motions of wavelength $\lambda _z = 2\pi /\kappa _z$ to the energy of the flow.

3.3. Triple-decomposed momentum and velocity budgets

Additionally to the spectra, the energy budget of the streamwise dispersive velocity component $\tilde {u}$ will be analysed at a steady state in § 4. For completeness, and to shed light on the way energy is redistributed between the mean, dispersive and fluctuation fields, each of the corresponding budget equations will be presented in the following discussion (including equations that will not be further discussed in the paper). These budget equations can be easily obtained starting from the Reynolds-averaged momentum budget and from the budget equation of the Reynolds stress tensor (see, e.g., Davidson Reference Davidson2015):

(3.9) \begin{align} \frac {\partial \left \langle {u_i} \right \rangle }{\partial t} + \left \langle {u_k} \right \rangle \frac {\partial }{\partial x_k} \left \langle {u_i} \right \rangle & + \frac {1}{\rho }\frac {\partial \left \langle {P} \right \rangle }{\partial x_i} = \nu \nabla ^2 \left \langle {u_i} \right \rangle - \frac {\partial }{\partial x_k}\left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle .\end{align}

(3.10) \begin{align} \frac {\partial \big \langle {u^{\prime}_i u^{\prime}_j} \big \rangle }{\partial t} + \left \langle {u_k} \right \rangle \frac {\partial }{\partial x_k} \big \langle {u^{\prime}_i u^{\prime}_j} \big \rangle \,=\, & - \,\left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle \frac {\partial \langle {u_j} \rangle }{\partial x_k} - \big \langle {u^{\prime}_ju^{\prime}_k} \big \rangle \frac {\partial \left \langle {u_i} \right \rangle }{\partial x_k} \nonumber\\ & - \,\frac {\partial }{\partial x_k} \left \langle {u^{\prime}_i u^{\prime}_j u^{\prime}_k} \right \rangle - \frac {\partial }{\partial x_i} \left \langle {u^{\prime}_j\frac {P^{\prime}}{\rho }} \right \rangle - \frac {\partial }{\partial x_j} \left \langle {u^{\prime}_i \frac {P^{\prime}}{\rho }} \right \rangle \nonumber\\ & + \,\left \langle {\frac {P^{\prime}}{\rho }\left (\frac {\partial u^{\prime}_i }{\partial x_j} + \frac {\partial u^{\prime}_j}{\partial x_i}\right )} \right \rangle \nonumber \\ & + \,\nu \nabla ^2 \big \langle {u^{\prime}_i u^{\prime}_j} \big \rangle - 2\nu \left \langle {\frac {\partial u^{\prime}_i }{\partial x_k}\frac {\partial u^{\prime}_j}{\partial x_k}} \right \rangle .\end{align}

The above equations are usually derived by time-averaging, whereas the present work resorts to averaging in the streamwise direction, over multiple repetitions of a simulation and over multiple phases of a periodic domain. Nevertheless, it is trivially shown that the above equations are valid regardless of how averaging is performed as long as the averaging operator fulfils the following properties (as it does in this case):

(3.11) \begin{equation} \left \langle {f^\prime } \right \rangle = 0; \,\,\, \left \langle {\left \langle {f} \right \rangle } \right \rangle = \left \langle {f} \right \rangle ; \,\,\, \left \langle {\left \langle {f} \right \rangle f^\prime } \right \rangle = \left \langle {f} \right \rangle \left \langle {f^\prime } \right \rangle ,\end{equation}

where $f$ indicates a generic random function. Equation (3.9) is only valid for the expected velocity field $ \langle {\boldsymbol {u}} \rangle$ ; leveraging the fact that the triple decomposition is a particular case of the Reynolds one, a budget equation for the $i$ th component $U_i$ of the mean field can be obtained by substituting (3.1) in (3.9) and taking a spanwise average:

(3.12) \begin{equation} \frac {\partial U_i}{\partial t} + U_k\frac {\partial U_i}{\partial x_k} + \frac {1}{\rho }\frac {\partial \left \langle {P} \right \rangle _z}{\partial x_i} = \nu \nabla ^2 U_i - \frac {\partial }{\partial x_k}\left \langle {\tilde {u}_i\tilde {u}_k} \right \rangle _z - \frac {\partial }{\partial x_k} \left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle _z. \end{equation}

The above equation resembles the Reynolds-averaged Navier–Stokes momentum equation (3.9), except that an additional term appears. Not only does the mean field feel the presence of turbulence through the Reynolds stress $ \langle {u^{\prime}_i u^{\prime}_k} \rangle _z$ , but it is also influenced by the dispersive field through the dispersive stress $ \langle {\tilde {u}_i\tilde {u}_k} \rangle _z$ . Notice, moreover, that (2.4), which describes the time-evolution of the bulk velocity, can be obtained by integrating (3.12) over the flow domain. By once again substituting (3.1) in (3.9) and by subtracting (3.12), one obtains a balance equation for the dispersive momentum:

(3.13) \begin{align} \frac {\partial \tilde {u}_i}{\partial t} + \left (U_k + \tilde {u}_k\right )\frac {\partial \tilde {u}_i}{\partial x_k} + \frac {1}{\rho }\frac {\partial \tilde {p}}{\partial x_i} = \nu \nabla ^2\tilde {u}_i &\underbrace {- \tilde {u}_k\frac {\partial U_i}{\partial x_k}}_{\text {V4}} + \frac {\partial }{\partial x_k} \left \langle {\tilde {u}_i \tilde {u}_k} \right \rangle _z \\ &- \frac {\partial }{\partial x_k} \left ( \left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle - \left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle _z \right ) \notag, \end{align}

where $\tilde {u}_i$ indicates the $i$ th component of $\boldsymbol {u}_d$ . Notice that (3.12) and (3.13) are of general validity (symmetries and simplifications due to the geometry have not been considered yet). Starting from the momentum budgets (3.12) and (3.13), energy budgets are trivially obtained by multiplication with $U_i$ and $\tilde {u}_i$ , respectively, and by performing a subsequent spanwise average. Additionally, an energy budget for the fluctuation field can be obtained by substituting (3.1) in (3.10) and by spanwise-averaging. After rearranging some terms, one obtains

(3.14) \begin{align} \frac {\partial }{\partial t}\frac {U_i^2}{2} \,+\, U_k \frac {\partial }{\partial x_k}\frac {U_i^2}{2} \,=\,\, &\nu \nabla ^2\frac {U_i^2}{2} \,-\, \nu \left (\frac {\partial U_i}{\partial x_k}\frac {\partial U_i}{\partial x_k}\right ) \,-\, \frac {U_i}{\rho }\frac {\partial \left \langle {P} \right \rangle _z}{\partial x_i} \\ \notag & + \underbrace {\left \langle {\tilde {u}_i\tilde {u}_k} \right \rangle _z\frac {\partial U_i}{\partial x_k}}_{\text {M5}} \,+\, \left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle _z\frac {\partial U_i}{\partial x_k} \\ \notag &- \frac {\partial }{\partial x_k}\left [ U_i \left ( \left \langle {\tilde {u}_i\tilde {u}_k} \right \rangle _z + \left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle _z \right ) \right ], \end{align}

(3.15) \begin{align} \frac {\partial }{\partial t}\frac {\left \langle {\tilde {u}_i^2} \right \rangle _z}{2} \,+\, \underbrace {\frac {\partial }{\partial x_k}\left \langle { (U_k + \tilde {u}_k) \frac {\tilde {u}_i^2}{2} } \right \rangle _z}_{\text {D1}} \,=\,\, &\underbrace {\phantom {\Big |\!}\nu \nabla ^2\frac {\left \langle {\tilde {u}_i^2} \right \rangle _z}{2}}_{\text {D2}} \,\underbrace {-\, \nu \left \langle {\frac {\partial \tilde {u}_i}{\partial x_k}\frac {\partial \tilde {u}_i}{\partial x_k}} \right \rangle _z}_{\text {D3}} \,\underbrace {-\, \left \langle {\frac {\tilde {u}_i}{\rho }\frac {\partial \tilde {P}}{\partial x_i}} \right \rangle _z}_{\text {D4}} \\ \notag & \underbrace {- \left \langle {\tilde {u}_i \tilde {u}_k} \right \rangle _z \frac {\partial U_i}{\partial x_k}}_{\text {D5}} \,+\, \underbrace {\left \langle {\left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle \frac {\partial \tilde {u}_i}{\partial x_k}} \right \rangle _z}_{\text {D6}}\\ \notag & \underbrace {- \frac {\partial }{\partial x_k} \left \langle {\tilde {u}_i\left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle } \right \rangle _z}_{\text {D7}}, \end{align}

(3.16) \begin{align} \frac {\partial }{\partial t} \frac {\left \langle {u^{\prime}_i u^{\prime}_i } \right \rangle _z}{2} + \frac {\partial }{\partial x_k}\!\left \langle {\!\left (U_k + \tilde {u}_k\right )\! \frac {\left \langle {u^{\prime}_i u^{\prime}_i } \right \rangle }{2}\!} \right \rangle _z\! =\,\, & + \,\nu \nabla ^2 \frac {\left \langle {u^{\prime}_i u^{\prime}_i } \right \rangle _z}{2} \,-\, \nu \left \langle {\frac {\partial u^{\prime}_i }{\partial x_k}\frac {\partial u^{\prime}_i }{\partial x_k}} \right \rangle _z \notag\\ & - \left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle _z\frac {\partial U_i}{\partial x_k} \,\underbrace {-\, \left \langle {\left \langle {u^{\prime}_i u^{\prime}_k} \right \rangle \frac {\partial \tilde {u}_i}{\partial x_k}} \right \rangle _z}_{\text {F6}} \notag\\ & -\frac {1}{2}\frac {\partial }{\partial x_k} \left \langle {u^{\prime}_i u^{\prime}_i u^{\prime}_k} \right \rangle _z \,-\, \frac {\partial }{\partial x_i} \left \langle {\frac {u^{\prime}_i }{\rho }P^{\prime}} \right \rangle _z. \end{align}

The above equations were also found by Reynolds & Hussain (Reference Reynolds and Hussain1972) using a slightly different averaging technique. Consider the equation for the dispersive kinetic energy (3.15). The dispersive kinetic energy gets transported by both the mean and the dispersive fields (term D1); a pressure term appears in the equation (D4), as well as the usual viscous diffusion (D2) and dissipation (D3) terms. The former indicates that viscosity tends to smear the dispersive energy out over time, whereas the latter represents the power lost to viscous forces. The two remaining equations ((3.14) and (3.16)) all share analogous transport, pressure and viscous terms; the only difference is that the mean kinetic energy $U_i^2$ is only transported by the mean field $U_i$ , and not by the dispersive one $\tilde {u}_i$ .

Most importantly, the above equations shed a light on how the mean, dispersive and fluctuation fields exchange energy. For instance, the dispersive stress term D5 appears both in the balance of dispersive energy (3.15) and in the balance of mean energy (3.14) (term M5) with opposite sign: it thus represents an exchange of power between the mean and the dispersive fields. Similarly, term D6 represents an exchange of power between the dispersive and fluctuation fields, as it appears both in the dispersive balance and in the fluctuation one (3.16) (term F6) with opposite sign. To sum up, the dispersive stresses enable the exchange of energy between the mean and the dispersive fields, whereas the Reynolds stresses allow the exchange of energy between the dispersive and the fluctuation field (and additionally between the mean and the fluctuation field, as is usual; see (3.14) and (3.16)).

In § 4, only the energy budget of the streamwise dispersive energy $\tilde {u}^2$ will be analysed; indeed, the quantities involved in the budgets of the two remaining components $\tilde {v}^2$ and $\tilde {w}^2$ are too small compared with the fluctuations to be captured with a satisfactory signal-to-noise ratio. Such a $\tilde {u}^2$ -budget can be obtained in a similar way to (3.15), but without averaging in the spanwise direction; after considering all the simplifications and symmetries due to the present geometry, the budget reads

(3.17) \begin{align} \frac {\partial }{\partial t} \frac {\tilde {u}^2}{2} + \underbrace {\left ( \tilde {v}\frac {\partial }{\partial y} + \tilde {w}\frac {\partial }{\partial \zeta } \right ) \frac {\tilde {u}^2}{2}}_{\mathcal {T}_a} \,=\, &\underbrace {\nu \nabla ^2\frac {\tilde {u}^2}{2} -\nu \left (\frac {\partial \tilde {u}}{\partial y}\frac {\partial \tilde {u}}{\partial y} + \frac {\partial \tilde {u}}{\partial \zeta }\frac {\partial \tilde {u}}{\partial \zeta } \right )}_{\mathcal {V}} + \underbrace {\tilde {u}\frac {\partial }{\partial y}\left \langle {\tilde {u}\tilde {v}} \right \rangle _z}_{\mathcal {T}_c} \nonumber\\ &\underbrace {-\,\tilde {u}\tilde {v}\frac {\partial U}{\partial y}}_{\mathcal {P}} \underbrace {-\, \tilde {u}\frac {\partial }{\partial y}\left (\left \langle {u^{\prime}v^{\prime}} \right \rangle - \left \langle {u^{\prime}v^{\prime}} \right \rangle _z \right )}_{\mathcal {T}_{uv}} \underbrace {- \tilde {u}\frac {\partial }{\partial \zeta }\left \langle {u^{\prime}w^{\prime}} \right \rangle }_{\mathcal {T}_{uw}}. \end{align}

An additional term appears here with respect to (3.15) – that is, $\mathcal {T}_c$ . This term integrates to zero when spanwise-averaged, both explaining its absence from (3.15) and indicating that the term only spatially redistributes $\tilde {u}$ -energy. The remaining terms all have a correspondent in (3.15). Notice that here, both the viscous diffusion (D2 in (3.15)) and the viscous dissipation (D3) are grouped in a single term $\mathcal {V}$ . The terms $\mathcal {T}_{uv}$ and $\mathcal {T}_{uw}$ quantify the work done by the Reynolds stresses on the dispersive velocity $\tilde {u}$ . Such work was rewritten in (3.15) as two separate terms (D6 and D7), of which one (D6) corresponds to a power exchange with the fluctuation field, and the other (D7) can be easily shown to turn zero for the $\tilde {u}$ component when integrated over the flow domain (notice that the same does not hold for the $w$ component owing to the slip length boundary condition). That is, term D7 only yields a spatial redistribution of $\tilde {u}$ -energy.

Finally, the $\mathcal {P}$ term of (3.17) is of particular interest for the results of this paper. We refer to such a term as dispersive production owing to its similarity to the canonical turbulence production term in a channel flow (Davidson Reference Davidson2015). Such a term corresponds to the D5 term in (3.15), and it thus represents an energy exchange with the mean field. The term accounts for the work done by the momentum flux V4 in (3.13) on the velocity $\tilde {u}$ . The momentum flux is better discussed by considering the balance of streamwise momentum and by applying all simplifications due to the geometry:

(3.18) \begin{equation} -\tilde {v}\frac {\partial U}{\partial y}. \end{equation}

It is clear from (3.18) that the discussed momentum flux originates from $\tilde {v}$ transporting the mean field $U$ . In other words, the circulatory dispersive $\tilde {v}$ - $\tilde {w}$ motions do not directly provide energy to the streamwise component $\tilde {u}$ , but still passively enable the transfer of energy from the mean flow to $\tilde {u}$ by transporting $U$ -momentum. A more intuitive explanation of how the $\tilde {u}$ field is indirectly produced by the circulatory motions (under certain circumstances) will be provided in § 4 using flow visualisations as an example.

3.4. Time scale for the decay of secondary motions; volume and plane averages

In § 5, the time needed by secondary motions to decay will be measured as per the objective of this study. Before doing so, such a time scale needs to be defined. Defining a time scale for the decay of turbulent eddies is a largely subjective process. For instance, Flores & Jiménez (Reference Flores and Jiménez2010) found that the log-layer of turbulent flows in a restricted simulation box bursts quasiperiodically, and linked the estimated period to the life span of log-layer eddies. LeHew et al. (Reference LeHew, Guala and McKeon2013) and Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014), instead, resorted to identifying turbulent coherent structures and tracking them in time; their lifespan is given by the distance in time between their first and last identification. In our case, no sophisticated strategy is needed to track the secondary motions, as their spatial position is fixed and their features are satisfactorily captured by the dispersive velocity field $\boldsymbol {u}_d$ (as will be discussed in § 4). We thus define some energy measure $e(t)\geqslant 0$ using the dispersive velocity and track it in time; the energy will start from a value $e(0)$ seen at the steady state and then decline to zero for $t\to \infty$ as the secondary motions decay. The time scale $T$ for the decay can then be defined as the time required for most of the energy to vanish; more precisely, as the minimum value of $t$ after which the energy $e$ never exceeds a threshold $\epsilon$ :

(3.19) \begin{eqnarray} T \,\mid \, e(T)=\epsilon ,\,\, e(T+\Delta t) \lt \epsilon \,\,\,\forall \,\,\, \Delta t \gt 0. \end{eqnarray}

A graphical representation of the above definition is provided in figure 2; in the figure, the generic energy measure $e$ is replaced by the volume-averaged quantity $I_u$ (defined in the next paragraph). The threshold $\epsilon$ is set to $15\, \%$ of the initial value $e(0)$ for multiple reasons. First, having a large threshold is beneficial for the signal-to-noise ratio: the smaller the measured value of $e$ , the greater its relative statistical uncertainty. It is moreover desirable to let the value of the threshold depend on the initial condition: secondary motions with different values of $\Lambda _s/h$ hold different amounts of energy at a steady state (Wangsawijaya & Hutchins Reference Wangsawijaya and Hutchins2022). It would not be sensible, then, to compare the time to decay of secondary motions of different sizes as measured by a fixed $\epsilon$ : every secondary motion starts from a different initial energy value. Instead, by letting $\epsilon = 15\, \% \,e(0)$ , the inverse $1/T$ of the time scale measures some sort of generalised decay rate ( $1/T$ would exactly be a multiple of the decay rate in the case of an exponential decay).

It remains to be specified which energy measure is used to calculate the time scale. In this respect, two different approaches will be pursued. To begin with, the energy of the dispersive field will be volume-averaged to yield two values $I_u$ and $I_{vw}$ for the streamwise and circulatory patterns, respectively:

(3.20) \begin{align} I_{u}(t) & = \frac {1}{2h}\int _0^{2h} \frac {1}{2} \big\langle {\tilde {u}^2} \big\rangle _z \,\text {d}y\text {,}\end{align}

(3.21) \begin{align} I_{vw}(t) & = \frac {1}{2h}\int _0^{2h} \frac {1}{2} \big\langle {\tilde {v}^2 + \tilde {w}^2} \big\rangle _z \,\text {d}y\text {.} \end{align}

By applying the (3.19) to the above quantities, the time scales $T_u$ and $T_{vw}$ for the streamwise pattern and for the circulatory motions are obtained. The streamwise pattern and the circulatory motions are indeed two separate features of the dispersive fields as discussed in § 3.1; they will be thus treated separately.

Figure 2. Graphical representation of the definition of the time $T_u$ needed for the pattern of streamwise dispersive velocity to decay. Such a time scale is defined by applying (3.19) to the volume-averaged energy $I_u$ . The dark green line represents the time evolution of $I_u$ ; the lighter green area its $95\, \%$ confidence interval (as per § 3.5). Similarly, a marker with an error bar is used to indicate $T_u$ and its $95\, \%$ confidence interval (see, once again, § 3.5). The horizontal dotted line indicates the threshold value $\epsilon =0.15\, I_u(0)$ used to define $T_u$ . Data at $Re_\tau =500$ , $\Lambda _s/h = 2$ .

The above approach defines two scalar quantities – $I_u$ and $I_{vw}$ – which can be easily tracked in time; their interpretation is also straightforward, as they represent the cumulative amount of energy held by secondary motions. However, such an approach is unable to capture the spatial complexity of the secondary motions. To recover some of it, for instance, one might resort to averaging the dispersive energy on wall-parallel planes. As an example, the plane-average $i_u$ of the streamwise energy is defined as follows:

(3.22) \begin{equation} i_u(t,y) = \frac {1}{2}\big \langle {\tilde {u}^2} \big \rangle _z. \end{equation}

Doing so, wall-normal energy variations can be investigated. Bear in mind (see (3.8)) that the value of $i_u$ seen at some wall-normal position $y$ is equal to the sum of all the energy contributions (as measured by the spectrum $\Phi _{\tilde {u}\tilde {u}}$ ) of the Fourier modes $\mathcal {F}_z \{\tilde {u}\}$ at that same position $y$ . Leveraging this property, one can also apply (3.19) of the time scale to the dispersive spectrum $\Phi _{\tilde {u}\tilde {u}}$ . The resulting time scale describes the life time of motions of given spanwise wavelength $\lambda _z = 2\pi /\kappa _z$ at a given wall-normal height $y$ .

3.5. Estimation of the dispersive field and uncertainty quantification

As previously explained, expected values are indicated by $ \langle {\cdot } \rangle$ throughout the present article. With the exception of this section, the same symbol is also used to indicate the estimates of such expected values as computed using the present data. The estimates deviate from the actual expected values as they are calculated on a sample of finite size; in other words, the estimates are affected by statistical uncertainty. It is assumed that statistical uncertainty dominates other sources of error (such as discretisation and round-off error, or the error introduced by the limited domain size), so that the overall uncertainty on the estimated statistics is given by the statistical uncertainty alone. This section deals with the quantification of such statistical uncertainty.

In particular, the uncertainty on the dispersive velocity field is of interest – as its values are referenced in the main results of this study. For the sake of clarity, the estimator of the dispersive velocity will be indicated by $\breve {u}_i$ , whereas its exact (theoretical) value will be $\tilde {u}_i$ as usual:

(3.23) \begin{equation} \tilde {u}_i \approx \breve {u}_i. \end{equation}

Here, the subscript $i$ indicates the component, so that $(\tilde {u}_1, \tilde {u}_2, \tilde {u}_3) = (\tilde {u}, \tilde {v}, \tilde {w})$ . The dispersive field is estimated through a streamwise-, phase- and ensemble-average; such operation can be reinterpreted as the arithmetic mean of a set of intermediate averages $\overline {u}^{x\Delta }_{i,k}$ :

(3.24) \begin{equation} \breve {u}_i = \frac {1}{N_s} \sum _{k=1}^{N_s} \,\underbrace {\frac {\Lambda _s}{L_z}\int _{0}^{L_z} \Delta _z(z+\zeta )\, \left ( \frac {1}{L_x} \int _{0}^{L_x} u_{i,k}-\big \langle {u_{i,k}} \big \rangle _z\, \text {d}x\right )\, \text {d}z}_{\overline {u}^{x\Delta }_{i,k}}, \end{equation}

where the index $k$ refers to one of the $N_s$ repetitions of a given simulation and $\Delta _z$ is a Dirac comb function used for the computation of the phase-average; by indicating Dirac’s delta distribution as $\delta$ ,

(3.25) \begin{equation} \Delta _z = \sum _{j=1}^{L_z/\Lambda _s - 1} \delta \left (z+j\Lambda _s\right ). \end{equation}

The $N_s$ values of $\overline {u}^{x\Delta }_{i,k}$ computed from the many repetitions of each simulation should be reasonably uncorrelated owing to the discussion of § 2.1. In light of the central limit theorem (Billingsley Reference Billingsley1995), then, the probability distribution function of $\breve {u}_i$ can be modelled by a normal distribution (indicated by $\mathcal {N}$ ) with mean $\tilde {u}_i$ ; its standard deviation $\breve {\sigma }_i$ can be estimated from the standard deviation $\overline {\sigma }^{x\Delta }_i$ of the set of values of $\overline {u}^{x\Delta }_{i,k}$ (which, in turn, is directly computed):

(3.26) \begin{align} \breve {u}_i & \sim \mathcal {N}\left (\tilde {u}_i, \breve {\sigma }_i\right ), \end{align}

(3.27) \begin{align} \breve {\sigma }_i & = \frac {\overline {\sigma }^{x\Delta }_i}{\sqrt {N_s}}.\end{align}

Note that, in the context of the present discussion, the cross-correlation between different components of $\breve {u}_i$ is neglected; this simplification does not affect the estimates of the uncertainty on $I_u$ . Owing to the above discussion, the $95\, \%$ confidence interval for the estimate of $\tilde {u}_i$ is given by

(3.28) \begin{equation} \tilde {u}_i = \breve {u}_i \pm \mathcal {E}\!\left \{\tilde {u}_i\right \} = \breve {u}_i \pm 2\breve {\sigma }_i. \end{equation}

Hence, the uncertainty $\mathcal {E}\!\{\tilde {u}_i\}$ can be propagated, as an example, to the energy $\tilde {u}^2/2$ of the streamwise field:

(3.29) \begin{equation} \mathcal {E}\!\left \{\frac {\tilde {u}^2}{2}\right \} \,\approx \, \sqrt {\left ( \frac {\text {d} \tilde {u}^2/2}{\text {d} \tilde {u}} \right )_{\breve {u}}^2 \left (\mathcal {E}\!\left \{\tilde {u}\right \}\right )^2 } \,=\, |\breve {u}|\,\mathcal {E}\!\left \{\tilde {u}\right \}. \end{equation}

Next, the uncertainty is propagated to the estimate of $I_u$ . In theory, this would require modelling the correlation between the values of $\tilde {u}^2/2$ at different spatial position. However, a pessimistic estimate of the uncertainty on $I_u$ can be found as follows. It is assumed that extreme events (meaning events for which $(\overline {u}^{x\Delta })^2/2$ falls outside of the $95\, \%$ confidence interval for $\tilde {u}^2/2$ ) happen simultaneously at all spatial positions. In other words, given that an extreme event occurs at some spatial position, one is certain to observe an extreme event at any other spatial position. This is a pessimistic assumption: likely, extreme events are confined in space and do not involve the whole flow domain. However, such an assumption might account for large eddies possibly triggering extended coherent regions of extreme events. In light of the above assumption, the error on $I_u$ is estimated as

(3.30) \begin{equation} \mathcal {E}\!\left \{I_u\right \} \,\approx \, \frac {1}{2h} \int _0^{2h} \left \langle {\mathcal {E}\!\left \{\frac {\tilde {u}^2}{2}\right \}} \right \rangle _z \,\text {d}y. \end{equation}

Similar considerations can be leveraged to estimate the uncertainty affecting $I_{vw}$ and $\Phi _{\tilde {u}\tilde {u}}$ .

Finally, the uncertainty can be propagated to the time scale $T_u$ defined by $I_u$ . To do so, the (estimated) value of the time derivative $I_{u,t}$ of $I_u$ at the time $T_u$ is used:

(3.31) \begin{equation} \mathcal {E}\!\left \{T_u\right \} \approx \sqrt {\frac {1}{I_{u,t}^2(T_u)}}\, \mathcal {E}\!\left \{I_u\right \}. \end{equation}

Once again, a similar procedure can be used to quantify the uncertainty on the time scales defined by $I_{vw}$ and $\Phi _{\tilde {u}\tilde {u}}$ .

5.1. Volume-averaged dispersive energy

As discussed in section § 3.4, the time needed for the dispersive field (or, to be more precise, by some of its features) to decay can be quantified by defining some energy measure and by tracking it in time. When the rough strips are removed from the walls of the channel, the dispersive velocity field starts losing energy with respect to the steady-state condition over strip-type roughness. A time scale for the decay can be then defined as the time needed by a generic energy measure to lose $85\, \%$ of its initial value. In particular, the volume-averaged energy of the momentum pathways ( $I_u$ , see (3.20)) and of the circulatory motions ( $I_{vw}$ , see. (3.21)) will be used in this section to calculate two corresponding time scales ( $T_u$ and $T_{vw}$ , see § 3.4 and (3.19)). Using volume-averaged quantities yields easy-to-interpret results at the cost of neglecting the spatial complexity of the flow; an attempt to recover some of this spatial complexity will be later done in § 5.2 and in Appendix B.

The initial steady-state values of $I_u$ and $I_{vw}$ are shown in figure 8(a–d) for each available flow configuration; panels (e–h) show their time evolution. At the lower Reynolds number (figure 8 a,b), motions with $1 \leqslant \Lambda _s/h \leqslant 4$ hold a similar amount of energy independently of their size. This contradicts the observation (Wangsawijaya & Hutchins Reference Wangsawijaya and Hutchins2022) that structures of period $\Lambda _s \approx 2h$ ( $s\approx h$ ) are more energetic than those of any other size; we suggest this discrepancy to be a consequence of the low Reynolds number. Indeed, low- $Re$ data suffers from two issues. It is known, for instance, that linear transient growth analysis predicts large structures to show significant values of transient growth only for sufficiently high Reynolds numbers (Cossu et al. Reference Cossu, Pujals and Depardon2009); moreover, large-scaled outer-layer eddies only become significantly energetic in channel flows if the Reynolds number is high enough (Lee & Moser Reference Lee and Moser2015). It can be thus expected that if the flow acts to favour structures of a specific outer-scaled size, this would only be seen at a sufficiently high Reynolds number. Our data at $Re_\tau =500$ (panels c,d) confirm this line of reasoning: the energy held by secondary motions is maximum for $\Lambda _s = 2h$ as expected. A second issue with low- $Re$ data is given by the lack of scale separation: at the considered Reynolds number ( $Re_\tau = 180$ ), the outer-layer length scale $h$ is equivalent to $180\,\delta _{v}$ . Such a value is not too far from the dominant spanwise scale of near-wall structures ( $100\,\delta _{v}$ , see Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967). In other words, we cannot easily tell whether the secondary motions that we impose are inner- or outer-scaled. Care must be then exerted in interpreting low- $Re$ data.

Figure 8. Volume-averaged dispersive energy. Initial (steady-state) values of the volume-averaged energy for varying spanwise period $\Lambda _s$ of the roughness pattern: (a) $I_u^+$ , $Re_\tau =180$ ; (b) $I_{vw}^+$ , $Re_\tau =180$ ; (c,d) same as panels (a,b), but at $Re_\tau =500$ . The error bars indicate the $95\, \%$ confidence interval, estimated as per § 3.5. Time evolution of the volume-averaged energy normalised by its initial value: (e) $I_u(t)/I_u(0)$ , $Re_\tau =180$ ; (f) $I_{vw}(t)/I_{vw}(0)$ , $Re_\tau =180$ ; (g–h) same as panels (e–f), but at $Re_\tau =500$ . The horizontal dotted line indicates the threshold value used for the calculation of the time scale (see (3.19)). Colour legend as in table 1: $\Lambda _s/h = 0.5$ ; $\Lambda _s/h = 1$ ; $\Lambda _s/h = 2$ ; $\Lambda _s/h = 4$ ; $\Lambda _s/h = 6$ .

As for the time evolution of the volume averages (figure 8 e–h), they typically show a monotonically decreasing trend both for the streamwise and circulatory patterns. There are, however, notable exceptions. At low Reynolds number (panel e), the energy of the streamwise pattern temporarily exceeds its initial value for $\Lambda _s/h = 4$ , to then decay as expected; similarly, it temporarily increases after an initial decay for $\Lambda _s/h=6$ . Most importantly, at the higher Reynolds number (panel g) the volume-averaged energy monotonically decays in all cases, except for the energy of the streamwise pattern for $\Lambda _s/h=2$ and $\Lambda _s/h=4$ . An excess of streamwise energy with respect to the initial condition is seen for these cases. To better illustrate our findings, we define a transient growth coefficient $G_u$ :

(5.1) \begin{equation} G_u = max _t \frac {I_u(t)}{I_u(0)} \end{equation}

and plot it against the size of the secondary motions in figure 9(a,b). At $Re_\tau =180$ (panel a), transient growth is only observed for $\Lambda _s/h=4$ and is not particularly pronounced ( $3\, \%$ ). At $Re_\tau =500$ (Panel b), a modest ( $6\, \%$ ) transient growth of the dispersive $\tilde {u}$ -energy is seen for $\Lambda _s/h = 2$ ; a larger growth ( $13\, \%$ ) is seen for $\Lambda _s/h = 4$ . The confidence level on the occurrence of transient growth given the statistical uncertainty is estimated in Appendix A to be $78.9\, \%$ and $85.2\, \%$ for $\Lambda _s/h = 2$ and $\Lambda _s/h = 4$ , respectively, at $Re_\tau =500$ . No excess energy is seen for other values of $\Lambda _s$ at the higher Reynolds number, except for a negligible transient growth for $\Lambda _s/h = 1$ . It appears that, for $2 \leqslant \Lambda _s \leqslant 4$ , the generation of the high- and low-momentum pathways continues for a short time even after removing the rough strips that enable their sustainment. This might be evidence in favour of the hypothesis of Townsend (Reference Townsend1976), who predicted that structures of this size ( $\Lambda _s \leqslant 4$ ) would be able to self-sustain: not only do the momentum pathways maintain their energy for a short time, but they even show excess energy with respect to the initial conditions. The excess energy is also reminiscent of linear transient growth analysis and its results (Del Álamo & Jiménez Reference Del Álamo and Jiménez2006), in the sense that we observe transient growth of $I_u$ before its decay. Moreover, the spanwise periods at which we observe transient growth ( $\Lambda _s/h \approx 2{-}4$ ) are in good agreement with the spanwise wavelength of maximum linear transient growth ( $\lambda _z/h=3$ ) found by the aforementioned study. There is, however, a crucial difference between our simulations and linear transient growth analysis, apart from the obvious nonlinearity of our system. Linear transient growth analysis studies the evolution of a perturbation of a given size, whereas we study the evolution of a fully developed structure.

Figure 9. Transient growth $G_u$ of the streamwise volume-averaged dispersive energy $I_u(t)$ ; (a) $Re_\tau =180$ , (b) $Re_\tau =500$ . The confidence level on the occurrence of transient growth is estimated in Appendix A. Time to decay (as defined per (3.19) applied to $I_u$ and $I_{uv}$ ) for the streamwise ( $T_u$ , black) and circulatory ( $T_{vw}$ , grey) dispersive energy; (c) $Re_\tau = 180$ , (d) $Re_\tau = 500$ . The dotted line represents a linear fit for $T_u$ performed (c) by rejecting data at $\Lambda _s/h=1$ and (d) by only considering data for $\Lambda _s/h\leqslant 2$ . Error bars indicate the $95\, \%$ confidence interval estimated as per § 3.5.

In linear transient growth analysis, the transient growth of $u$ -energy is found to be driven by the $v$ component through transport of the mean field (Del Álamo & Jiménez Reference Del Álamo and Jiménez2006). Similarly, in §§ 3.3 and 4, we have proposed that the formation of the $\tilde {u}$ momentum pathways might be driven by $\tilde {v}$ transport of the $U$ -field (dispersive production). It could be expected, then, that this same mechanism would drive the initial overshoot of $I_u$ energy that we observe. If this were the case, we would expect $I_{vw}$ to be able to maintain its initial energy for a longer time in cases for which excess $I_u$ is seen, or perhaps to decay slower than in other cases. A preliminary scrutiny of figure 8(f,h) suggests this is not the case: no initial plateau of $I_{vw}$ is seen for values of $\Lambda _s$ at which transient growth is observed. Moreover, we could not find any link between the rate of change of $I_{vw}$ at $t=0$ and the occurrence of transient growth. We further investigate the matter by calculating a time scale for the decay of $I_u$ and $I_{vw}$ based on the definition in (3.19). The resulting time scales are $T_{u}$ and $T_{vw}$ , respectively. Results are shown in figure 9(c,d) and are scaled with $h/u_p$ , as this time scale remains constant in physical terms across different simulations at the same $Re_\tau$ . Notice that one data point (the value of $T_u$ at $Re_\tau =500$ , $\Lambda _s/h=6$ , panel d) is affected by a significant amount of statistical uncertainty; the reasons for this can be multiple. On the one hand, this flow case might exhibit unexpectedly large velocity fluctuations, leading to a large statistical uncertainty on the dispersive velocity field (the uncertainty then propagates to $I_u$ and $T_u$ ). On the other hand, the dispersive velocity field is computed by performing, among others, a phase average over multiple spanwise periods (see § 2). For such a value of $\Lambda _s/h=6$ , the simulation box ( $L_z/h = 6$ ) only contains one period $\Lambda _s$ , so that no phase average is effectively performed. This might explain the increased statistical uncertainty of the flow case under consideration; in all other flow cases at the same Reynolds number, the simulation box contains at least two spanwise periods $\Lambda _s$ (see table 1) which can be exploited by the phase average. Finally, notice that the low- $Re$ simulation with matching $\Lambda _s$ (namely, $Re_\tau =180$ , $\Lambda _s/h=6$ , panel a) also has a relatively large statistical uncertainty; in this last case, nevertheless, it was possible to produce a larger number of repetitions of the simulation (as low- $Re$ simulations are cheap) so as to mitigate the statistical variability.

At the lower Reynolds number, $T_{vw}$ generally increases with $\Lambda _s$ , though the circulatory motions remain in the flow for an unexpectedly long time for $\Lambda _s/h = 1$ . At the higher Reynolds number, $T_{vw}$ increases with $\Lambda _s$ until it reaches a local maximum value for $\Lambda _s/h = 2$ . For larger values of $\Lambda _s$ , the circulatory motions remain confined in a $2h$ -wide region around roughness transitions (see figure 4), so that effectively they stop growing in size. Similarly, their time to decay remains bounded to $T_{vw}\leqslant 2$ (roughly). The local maximum of $T_{vw}$ for $\Lambda _s/h = 2$ might suggest that the longer-living circulatory motion is able to leverage production of $I_u$ for a longer time, causing the observed overshoot of $I_u$ . However, data at $\Lambda _s/h = 4$ contradict this idea: for this value of $\Lambda _s$ , $I_u$ undergoes an even stronger transient growth – and yet the circulatory motion decays in an unexpectedly short time. More generally, momentum pathways tend to live longer than the circulatory motions (as $T_u \gt T_{vw}$ in most cases). In §§ 3.3 and 4, we have argued that the momentum pathways are mostly produced by the circulatory motions through the dispersive production term: the observed delay in the decay of $I_u$ with respect to $I_{vw}$ supports this idea. Indeed, the streamwise dispersive field $\tilde {u}$ has its own dynamics, as indicated by (3.17); it is reasonable that it would be able to survive for a limited time after the circulatory motions have decayed. The absence of circulatory motions simply means that the $\tilde {u}$ field is not being fed energy by the dispersive production term anymore; the viscous dissipation will then gradually erode the remaining $\tilde {u}$ -energy until the decay is complete. In spite of the relevance of the dispersive production term, the variability of $T_{vw}$ fails to predict the variability of $T_u$ , indicating that production through $\tilde {v}$ -transport of $U$ alone cannot explain the dynamics of $I_u$ . Indeed, as once again highlighted in (3.17), the dynamics of $\tilde {u}$ is affected by several other terms.

The time to decay of $I_u$ generally linearly increases with the spanwise periodicity $\Lambda _s$ at the lower Reynolds number; at the higher $Re$ , it also generally increases, although the slope of the trend diminishes for $\Lambda _s/h \gt 2$ . There might be an exception to this monotonically increasing trend: the estimated value of $T_u$ for $\Lambda _s/h=6$ , $Re_\tau =500$ (figure 9 d) is essentially equal to that seen for $\Lambda _s/h=4$ at the same Reynolds number. However, data at $\Lambda _s/h=6$ are affected by a significant statistical uncertainty, making it impossible to assess whether the trend of $T_u$ effectively saturates, starts decreasing or keeps increasing (albeit with a smaller slope) for high $\Lambda _s$ . Interestingly, the observed transient growth of $I_u$ has no strong influence on $T_u$ : at $Re_\tau =500$ (panel d), the value of $T_u$ for $\Lambda _s/h = 2$ is either lower or comparable (accounting for the statistical uncertainty) than that for $\Lambda _s/h=6$ . Notice that transient growth is observed for $\Lambda _s/h = 2$ , but not for $\Lambda _s/h = 6$ . An additional noteworthy feature can be found in the low- $Re$ data of panel (c): momentum pathways of period $\Lambda _s = h$ subsist in the flow longer than what could be expected by extrapolating the observed trend. Although this could be a sign of the flow favouring these specific $h$ -scaled structures, it should be kept in mind that the structures might as well be inner-scaled owing to the lack of scale separation. As previously explained, a $1h$ period is indeed equivalent to $180\delta _{v}$ in wall units at $Re_\tau =180$ ; such an inner-scaled value of the period is not far from the typical spanwise spacing of buffer-layer small scales ( $100\delta _{v}$ , see Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967).

Starting from the temporal information we gather, we now estimate the streamwise distance needed for the $\tilde {u}$ pattern (that is, for the high- and low-momentum pathways) to decay. This enables the comparison of our results with other studies measuring the streamwise coherence of momentum pathways (e.g. Womack et al. Reference Womack, Volino, Meneveau and Schultz2022). First, we find a linear fit to the trend of $T_u$ against $\Lambda _s$ ; at $Re_\tau =180$ , we exclude the datapoint at $\Lambda _s/h=1$ , whereas at $Re_\tau =500$ , only the values $\Lambda _s/h \leqslant 2$ are considered. Although this is a rough approximation, the fit (shown by a dotted line in figure 9 c,d) correctly captures the order of magnitude of $T_u$ over the considered values of $\Lambda _s$ . Then, we estimate the streamwise distance $\Delta x_d$ needed for the decay of momentum pathways as $\Delta x_d = U_{b,s} \,T_u$ , where $U_{b,s}$ is the bulk velocity seen at a matching value of $Re_\tau$ between smooth walls. The underlying idea is that structures are advected downstream as they evolve and that the bulk velocity is the mean velocity at which this happens. Using the value of the bulk velocity seen between smooth walls is an arbitrary choice; again, it is only meant to roughly capture the order of magnitude of $\Delta x_d$ . The following crude estimates are obtained:

(5.2) \begin{align} \Delta x_d \approx 15.8 \,\Lambda _s \,\,\,\,\, \,\,\,\,\, &\text {at } Re_\tau = 180, \notag \\ \Delta x_d \approx 34.9 \,\Lambda _s - 12.7 h \,\,\,\,\, \,\,\,\,\, &\text {at } Re_\tau = 500\text{ for }\Lambda _s/h \leqslant 2. \end{align}

The order of magnitude of our estimates for a spanwise wavelength $\Lambda _s=h$ is $\Delta x_d \approx 15 - 20$ ; this is in line with the experiments of Womack et al. (Reference Womack, Volino, Meneveau and Schultz2022), who reported the streamwise extent of momentum pathways to be $18h$ at least.

As already discussed, the value of the time scale $T_{vw}$ associated with the circulatory motions saturates for high values of $\Lambda _s$ . One way of explaining the observed saturation could be given by the findings of Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014). The authors identified and tracked naturally occurring structures that are both wall-attached and associated with $u^{\prime}v^{\prime}$ -anticorrelation events; they found the lifetime $T_{ae}$ of these structures to be proportional to their half-height $h_y$ ,

(5.3) \begin{equation} T_{ae}^+ \approx 2h_y^+\text {.} \end{equation}

In the case of the present study, the time scale $T_u$ cannot be considered to be analogous to $T_{ae}$ , as it measures the lifespan of a velocity pattern that is not necessarily wall-attached and that includes regions in which $\tilde {u}$ and $\tilde {v}$ are positively correlated. The circulatory motions, instead, satisfy these two requirements: they extend upwards from the near-wall region and they are responsible for the $\tilde {u}$ - $\tilde {v}$ anticorrelation (see § 4). Thus, although the secondary motions are not features of the fluctuation field (differently from the structures observed by Lozano-Durán & Jiménez Reference Lozano-Durán and Jiménez2014), we will try to draw an analogy between $T_{vw}$ and $T_{ae}$ . By making (5.3) dimensional, we expect:

(5.4) \begin{equation} T_{vw} \approx 2 \frac {h_y}{u_p} \,\,\,\,\,\, \text {(expected)}. \end{equation}

To enable a comparison between our results and those of Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014), we now need to estimate the half-height $h_y$ of the circulatory motions. We define a vortex centre as the point at which, simultaneously, $\tilde {v}$ changes in sign in $z$ and $\tilde {w}$ changes in sign in $y$ . The so-defined vortex centres are marked in figure 4 and satisfactorily represent the centre of the cross-sectional circulatory motions. The structure half-height $h_y$ is then simply given by the wall-normal position of the vortex centre.

Figure 10. (a) Wall-normal position $h_y$ of the vortex centre of the circulatory motions against the spanwise period $\Lambda _s$ . (b) Time to decay $T_{vw}$ of the circulatory motions against the wall-normal position $h_y$ of the vortex centre. Light grey crosses indicate data at $Re_\tau = 180$ ; dark grey asterisks indicate data at $Re_\tau = 500$ . In panel (b), a linear fit to all available data is shown as a dashed line; error bars indicate the $95\, \%$ confidence interval, estimated as per § 3.5.

Figure 10(a) shows how the half-height $h_y$ changes with the spanwise period $\Lambda _s$ of the roughness pattern. A good collapse of low- and high- $Re$ data is observed. The circulatory motions get taller for increasing $\Lambda _s$ ; as $\Lambda _s$ approaches large values (say, $\Lambda _s/h \geqslant 2$ ), the growth of $h_y$ slows down. This saturation effect is expected: as can be seen from figure 4, the circulatory motions fill the entire channel half-height for $\Lambda _s/h\approx 2$ , so that their wall-normal growth is phyisically limited for larger values of $\Lambda _s$ . As the structures stop getting taller, they also become confined to a region surrounding roughness transitions (see § 4): it appears that the circulatory motions roughly maintain their $y$ - $z$ aspect ratio – in a way that is reminiscent of the attached eddy hypothesis (Marusic & Monty Reference Marusic and Monty2019). Finally, we compare our data to those of Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014) by plotting the time scale $T_{vw}$ of the circulatory motions against their half-height $h_y$ in figure 10(b): owing to (5.4), a linear trend is expected. Although the data do not appear to collapse on a single line, their dispersion appears to be directional, so that a linear fit is reasonable. By considering both low- and high- $Re$ data, we obtain

(5.5) \begin{equation} \frac {T_{vw,fit}}{h/u_p} \,=\, 4.48 \frac {h_y}{h} \,+\, 0.011. \end{equation}

Although the intercept of our fit is small (and thus in agreement with (5.3)), the slope we obtain is roughly two times larger than expected. This quantitative difference can be explained by the fact that we track different structures with respect to Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014); moreover, we use different definitions of the time scale and of the structure half-height – all of which are arbitrary. To conclude, the ansatz of (5.3) does not fully explain the variability of $T_{vw}$ , but it correctly captures the general trend of larger values of $h_y$ being typically associated with larger values of $T_{vw}$ . Then, the fact that the circulatory motions stop growing in size for large values of $\Lambda _s$ can help explain the observed saturation of $T_{vw}$ .

5.2. Plane-averaged dispersive energy: the decay of the wall shear stress pattern

Albeit easily interpretable, volume-averaged quantities hide the spatial complexity of the observed phenomena. To recover information in the wall-normal direction, we inspect the average $i_u$ (defined in (3.22)) of the streamwise dispersive energy on wall-parallel planes. Bear in mind that, as discussed in § 3.2, the value of $i_u$ seen at a given wall-normal position is equal to the sum of the contributions of each Fourier mode of the spectrum $\Phi _{\tilde {u}\tilde {u}}$ at the same distance from the wall.

Figure 11. Time evolution of the plane-averaged dispersive energy $i_u^+$ of the streamwise dispersive velocity. Please notice that, in spite of the logarithmic scale on the vertical axis, no premultiplication is used; consequently, the visual representation given by this figure is not well representative of the total (integral) amount of energy contained in different regions of the channel. All data at $Re_\tau =500$ ; (a) $\Lambda _s/h = 1$ , (b) $\Lambda _s/h=2$ and (c) $\Lambda _s/h=6$ .

The plane-averaged energy $i_u$ is shown in figure 11 for a selection of flow cases at $Re_\tau =500$ for a short time interval after the spanwise heterogeneity has been removed. Notice that no premultiplication is performed for this figure despite the log-scaled vertical axis: this is done to better highlight its near-wall features. As a consequence, though, the total amount of energy contained in the outer layer is misrepresented by the visualisation of figure 11. Similarly to what has been observed in figures 4 and 6, two separate energy peaks can be identified in figure 11 at the initial condition ( $t=0$ , steady state) for large enough values of $\Lambda _s/h$ : one in the viscous sublayer ( $y^+\leqslant 10$ ), one further away from the wall. This second peak is particularly pronounced for $\Lambda _s/h=2$ ( $y^+ \approx 200$ , panel b), and not particularly so for $\Lambda _s/h = 6$ ( $y^+ \approx 50$ , panel c); it is not observed at the lower Reynolds number or, at least, it is not as pronounced. In general, the streamwise pattern for $\Lambda _s/h = 6$ is much less energetic than that for $\Lambda _s/h=2$ , as previously observed by analysing volume averages (figure 8 c). The common feature of all data in figure 11 (which is also observed for values of $\Lambda _s$ and $Re_\tau$ that are not shown) is that the near-wall peak diffuses towards the core of the channel. Our steady-state analysis of the streamwise momentum budget (§ 4) suggests that viscous diffusion is responsible for this process. The near-wall peak is associated with the near-wall square-wave distribution of $\tilde {u}$ seen in figure 4, and hence to the spanwise distribution of wall shear stress.

Figure 12. Time to decay $T_w(\lambda _z)$ of the Fourier harmonics (represented by their wavelength $\lambda _z$ ) that compose the dispersive wall shear stress distribution, calculated by applying the definition in (3.19) to the spectrum $\Phi _{\tilde {u}\tilde {u}}$ (shown in figure 6) at $y^+ \approx 1$ . Only the first three odd energy-containing harmonics are shown. (a) $Re_\tau =180$ ; (b) $Re_\tau =500$ . The error bars indicate the $95\, \%$ confidence interval, estimated as per § 3.5. Colour legend as in table 1: $\Lambda _s/h = 0.5$ ; $\Lambda _s/h = 1$ ; $\Lambda _s/h = 2$ ; $\Lambda _s/h = 4$ ; $\Lambda _s/h = 6$ .

While the near-wall pattern of $\tilde {u}$ is particularly intense, it resides in a restricted region of the channel that shrinks in size as $Re_\tau$ increases. As a consequence, its contribution to the volume-averaged energy is marginal, so that the time scales measured in figure 9 are not representative of the decay of the dispersive wall shear stress distribution. A suitable time scale for this purpose could be defined, for instance, by applying the definition in (3.19) to the plane-averaged energy $i_u$ at the wall. Instead, we go one step further and analyse the individual contributions to $i_u$ of each Fourier mode (as per (3.8)) and measure their time to decay. In other words, we apply the definition in (3.19) of the time scale to the spectrum $\Phi _{\tilde {u}\tilde {u}}$ – evaluated at the first wall-normal grid point ( $y^+ \approx 1$ , please refer to table 1) to be representative of the wall shear stress. In this way, we obtain a different time scale $T_w(\lambda _z)$ for each Fourier mode (of wavelength $\lambda _z = 2 \pi /\kappa _z$ ) that constitues the pattern of wall shear stress.

Results are shown in figure 12, which collates data for all available values of $\Lambda _s/h$ . Only the first three energy-containing harmonics are shown; higher harmonics hold little energy so that the signal-to-noise ratio is excessively low. Patterns of wall shear stress of different period $\Lambda _s$ may contain Fourier modes with the same (or similar) spanwise wavelength $\lambda _z$ . The time needed for Fourier modes of comparable wavelengths from different simulations (different $\Lambda _s/h$ ) to decay is similar: data from different simulations at the same $Re_\tau$ appear to collapse on the same curve. In other words, it appears that the time evolution of the Fourier modes is influenced by their own wavelength $\lambda _z$ (and, of course, by the Reynolds number), but not much by the geometry of the problem (that is, by $\Lambda _s/h$ ) or by the amount of energy they hold at the initial steady state. As an example, the first harmonic of figure 6(a) and the third one in figure 6(d) hold a different amount of energy at the wall at the initial steady state; yet, they share the same spanwise wavelength $\lambda _z/h=2$ and take a comparable amount of time to decay as measured by $T_w$ (see figure 12).

Figure 12(a) shows data at $Re_\tau =180$ : the time to decay increases with the wavelength until $\lambda _z/h = 1$ to then saturate for larger values of $\lambda _z/h$ . At the higher Reynolds number (figure 12 b), instead, the time to decay increases up until its peak value at $\lambda _z/h = 2$ to then decrease for higher values of the wavelength. That is, Fourier modes of the wall shear stress with $\lambda _z/h = 2$ are the longest-lived ones. This value of the wavelength matches the period $\Lambda _s/h = 2$ for which secondary motions are the most energetic.

We interpret the above observations as follows. The longer time to decay might be a sign that near-wall Fourier modes of a specific wavelength ( $\lambda _z/h=2$ at $Re_\tau =500$ ) are less damped (e.g. by viscous dissipation) than the remaining modes during their evolution. This does not help to explain the generally increasing trend of $T_u$ against $\Lambda _s$ seen in figure 9(d): $T_u$ measures indeed the time to decay of a quantity that is integrated in the wall-normal direction over the entire channel half-height. Nevertheless, the near-wall behaviour discussed here might play a crucial role in the formation of secondary motions. Consider, for instance, a fully developed flow suddenly hitting a patch of spanwise heterogeneous roughness. The roughness can be idealised as a disturbance applied to the near-wall region; in light of the above discussion, we would expect disturbances of a specific size to be less damped than others in this region, so that they grow more energetic. This mechanism might help to explain why secondary motions with a spanwise period of $\Lambda _s/h = 2$ are the most energetic at $Re_\tau =500$ (see figure 8 c); bear in mind that a spanwise roughness pattern of period $\Lambda _s$ significantly excites the Fourier mode with wavelength $\lambda _z = \Lambda _s$ (as can be seen from the spectra of the roughness patterns in figure 6).