**2. Materials and Methods**

The present section is devoted to the methodology used for the numerical simulation of pressure-driven plane channel flow. The flow is governed by the incompressible Navier Stokes equations. Channel flow is described here using the Cartesian coordinates *x*,*y*,*z*, respectively the streamwise, wall-normal and spanwise coordinates. The velocity field *u*p*x*, *<sup>y</sup>*, *<sup>z</sup>*, *<sup>t</sup>*<sup>q</sup> is decomposed into the steady laminar base flow solution **<sup>U</sup>**p*y*q"p*Ux*, 0, 0<sup>q</sup> and a perturbation field *u*<sup>1</sup> p*x*, *y*, *z*, *t*q. Similarly, the pressure field is decomposed as *p*p*x*, *y*, *z*, *t*q " *xG* ` *p*<sup>1</sup> p*x*, *y*, *z*, *t*q. The equation governing the steady base flow for an incompressible fluid with constant density *ρ* and kinematic viscosity *ν* is given by

$$
\nu \frac{\partial^2 \mathcal{U}\_x}{\partial y^2} = \frac{1}{\rho} \mathcal{G} \tag{1}
$$

with *G* a constant. Together with the no-slip condition at the walls Equation (1) yields the analytic Poiseuille solution *Ux*91 ´ p*y*{*h*q2. The equation governing the perturbation field involves the base flow and reads

$$\frac{\partial \mathbf{u}'}{\partial t} + \mathbf{u}' \cdot \nabla \mathbf{u}' + \mathbf{U} \cdot \nabla \mathbf{u}' + \mathbf{u}' \cdot \nabla \mathbf{U} = -\frac{1}{\rho} \nabla p' + \nu \nabla^2 \mathbf{u}' \tag{2}$$

The channel geometry is formally infinitely extended, yet in the numerical representation it is given by its extent *Lx* ˆ 2*h* ˆ *Lz* as in Figure 1, with stationary walls at *y* " ˘*h* and periodic boundary conditions in *x* and *z*.

**Figure 1.** Schematic of the numerical domain with the laminar base flow profile (red).

The flow is driven by the imposed pressure gradient *G* assumed negative. The spanwise pressure gradient is explicitly constrained to be null. The centerline velocity *ucl* of the laminar base profile with the same pressure gradient is chosen as the velocity scale (*U*) and the half gap *h* of the channel is chosen as the lengthscale used for non-dimensionalization. Time is hence expressed in units of *h*{*U*. In these units the laminar velocity profile is given by *U*˚ *<sup>x</sup>* <sup>p</sup>*y*˚q " <sup>1</sup> ´ *<sup>y</sup>*<sup>2</sup> ˚. From Chapter 3 onwards only dimensionless quantities will be used and the ˚ notation will be dropped from there on. Primed quantities denote perturbations to the base flow while non-primed quantities involve the full velocity field, including the laminar base flow.

In the following we shall consider, both locally and temporally fluctuating quantities, as well as their time and space averages. We denote by ' the space (*x*, *<sup>z</sup>*) average and <sup>s</sup>'—the time average. Space-time averages are indicated by Ď¨. More explicitly the space-average operator is defined as the discrete average over the grid points, and the time average is the discrete average sum over the total number of snapshots in the steady regime.

Different velocity scales characterize the flow. One such scale is the centerline velocity *ucl* of the corresponding laminar flow with the same value of *G*. Another one is the total streamwise flow through the channel, *Ub* " Ě*ub*, where

$$u\_b(x, z, t) = \frac{2}{h} \int\_{-h}^{h} u\_x dy \tag{3}$$

is the so-called local bulk flow. Finally, the friction velocity is defined as *U<sup>τ</sup>* " p Ď*τ*{*ρ*q 1 <sup>2</sup> , where *τ* " p*τ<sup>t</sup>* ` *τb*q {2 ą 0, with *τ<sup>t</sup>* and *τ<sup>b</sup>* the net shear stress on the top and the bottom wall, respectively given by:

$$
\pi\_{t,b}(\mathbf{x}, z, t) = \pm \mu \frac{\partial \mathbf{u}\_x}{\partial y} \bigg|\_{t,b} \tag{4}
$$

where *μ* " *ρν* is the dynamic viscosity of the fluid. The three Reynolds numbers arising from these velocity scales are *Recl* " *uclh*{*ν*, *Reb* " *Ubh*{*ν* and *Re<sup>τ</sup>* " *Uτh*{*ν*. For the laminar base flow, they are inter-related as *Re*<sup>2</sup> *<sup>τ</sup>* " 3*Reb* " 2*Recl*. Imposing a pressure gradient *G* < 0 translates into a fixed average shear stress Ď*τ* on the walls which sets an imposed value of *Re<sup>τ</sup>* " *ReG <sup>τ</sup>* to stress that this is the control parameter.

Direct numerical simulation (DNS) of Equation (2) is carried out using the open source, parallel solver called Channelflow [38,39] written in C++. It is based on a Fourier–Chebychev discretization in space and a 3rd order semi-implicit backward difference scheme for timestepping. It makes use of the 2{3 dealiasing rule for the nonlinear terms. An influence matrix method is used to ensure the no-slip boundary condition at the walls. The numerical resolution is specified in terms of the spatial grid points p*Nx*, *Ny*, *Nz*q which translates into a maximum of p*Nx*{2 ` 1, *Nz*{2 ` 1q Fourier wavenumbers and *Ny* Chebychev modes. Please note that the definitions of *Nx* and *Nz* take into account the aliasing modes. The domain sizes used in this study, expressed in units of *h*, are *Lx* " 2*Lz* " 250 for 55 ă *ReG <sup>τ</sup>* ď 100 and *Lx* " 2*Lz* " 500 for 39 ď *ReG <sup>τ</sup>* ď 55. The local numerical resolution used is *Nx*{*Lx* " *Nz*{p2*Lz*q " 4.096 and *Ny* " 65, comparable to that used in Reference [34]. The simulation follows an "adiabatic descent": a first simulation is carried out at sufficiently high value of *ReG <sup>τ</sup>*, known to display space-filling turbulence. After the stationary turbulent regime is reached, *ReG <sup>τ</sup>* is lowered and the simulation advanced further in time. This step-by-step reduction has been performed down to *ReG <sup>τ</sup>* " 39. The initial condition for the first simulation is a random distribution of localized seeds of the kind described in Reference [40]. The time required *T* to reach a stationary regime gradually increases as *ReG <sup>τ</sup>* is decreased. As an order of magnitude, for *ReG <sup>τ</sup>* " 100, *T* « 1500, while for *ReG <sup>τ</sup>* " 50, *T* « 3000. Statistics are computed, after excluding such transients, from time series of lengths up to 2 ˆ 104 time units.
