# Input-output analysis of the stochastic Navier-Stokes equations: application to turbulent channel flow

Gilles Tissot, André V. G. Cavalieri, Etienne Mémin

► To cite this version:

Gilles Tissot, André V. G. Cavalieri, Etienne Mémin. Input-output analysis of the stochastic Navier-Stokes equations: application to turbulent channel flow. Physical Review Fluids, 2023, 8 (3), pp.1-18. 10.1103/PhysRevFluids.8.033904. [hal-04022896](https://HAL Id: hal-04022896)

Submitted on 10 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Input-output analysis of the stochastic Navier-Stokes equations:
application to turbulent channel flow

Gilles Tissot,¹, * André V. G. Cavalieri,² and Étienne Mémin¹

¹INRIA Rennes Bretagne Atlantique, IRMAR – UMR CNRS 6625, av. Général Leclerc, 35042 Rennes, France

²Department of Aerospace Engineering, Instituto Tecnológico de Aeronáutica,
Praça Mal. Eduardo Gomes 50, Vila das Acácias, 12228-900, São José dos Campos, Brazil

Stochastic linear modelling proposed in TISSOT, MÉMIN & CAVALIERI (J. Fluid Mech., vol. 912, 2021, A51) is based on classical conservation laws subject to a stochastic transport. Once linearised around the mean flow and expressed in the Fourier domain, the model has proven its efficiency to predict the structure of the streaks of streamwise velocity in turbulent channel flows. It has been in particular demonstrated that the stochastic transport by unresolved incoherent turbulence allows us to better reproduce the streaks through lift-up mechanism. In the present paper, we focus on the study of streamwise-elongated structures, energetic in the buffer and logarithmic layers. In the buffer layer, elongated streamwise vortices, named rolls, are seen to result from coherent wave-wave non-linear interactions, which have been neglected in the stochastic linear framework. We propose a way to account for the effect of these interactions in the stochastic model by introducing a stochastic forcing, which replaces the missing non-linear terms. In addition, we propose an iterative strategy in order to ensure that the stochastic noise is decorrelated from the solution, as prescribed by the modelling hypotheses. We explore the prediction abilities of this more complete model in the buffer and logarithmic layers of channel flows at $Re_τ = 180$, $Re_τ = 550$ and $Re_τ = 1000$. We show an improvement of predictions compared to resolvent analysis with eddy viscosity, especially in the logarithmic layer.

# I. INTRODUCTION

Coherent structures of the near-wall turbulence is an extensively explored topic. In the buffer layer, very close to the wall, the flow organises into streamwise vortices, or rolls, and elongated patterns of high/low streamwise velocity denoted as streaks [1–3]. These structures develop, break and are regenerated in a quasi-cyclical process [4]. A scenario explaining their behaviour [5] considers the cycle where the streaks intensify by the lift-up mechanism [6, 7], destabilise by spanwise meandering leading to non-linear interactions which give finally birth to new streamwise vortices. This final step allows to start a new cycle.

In the logarithmic layer, the flow organises as well along streaks [8, 9]. These structures are of larger size with a more disorganised motion due to the higher Reynolds number based on the wall distance (the wall distance in viscous units $y^+ = \frac{y u_τ}{\nu}$ is a Reynolds number based on the friction velocity $u_τ$, the kinematic viscosity $\nu$ and the wall distance in outer units $y$, typical length scale of the largest structure [10]). Understanding their dynamical behaviour is still today an active research area. Smaller scales appear to be unnecessary to sustain these structures [11], suggesting the presence of a self-sustaining mechanism at large scale. Several evidences indicate that this mechanism is similar to the one active in the buffer layer [12, 13]. As noted in Cossu and Hwang [12], these large scale coherent structures exist in the sense of (ensemble) averaging or filtering as associated to large eddy simulation (LES). Even if they do not need the small scales to survive, they still interact with them. As a consequence to predict their dynamical behaviour, it is crucial to include the effect of small scales on these large scales, through Reynolds stress models for instance, or, as we propose here, by stochastic modelling. As a practical example in Bae et al. [14], resolvent analysis has been used to extract these large coherent structures in view of performing diagnostics of their action in removing their contributions in a numerical simulation. This procedure yields a drastic reduction of the turbulence intensity. The reduction is significant in the buffer layer and slightly less so in the logarithmic-layer, highlighting the requirement of modelling improvements in this region. Besides, practical control strategies require an accurate prediction of these structures, and providing simplified models predicting coherent structures at a given scale with a high fidelity is still today challenging.

By the knowledge of the time-averaged velocity field and possibly of some higher-order statistics, predicting coherent structures in a turbulent flow without resolving the whole multiscale space-time dependent solution has become an important research direction, to which many groups have devoted strong efforts. Considering a linearisation of the Navier-Stokes operator around a suitably chosen flow – often taken as the time-averaged flow [15] – it is natural to

<footnote>* gilles.tissot@inria.fr</footnote>
search for wave-solutions in the Fourier domain, which beyond a natural physical meaning gives access to efficient
linear-algebra techniques. Since turbulence interacts with these wavy coherent structures, linearised solutions are
often insufficient, and a closure is required.

Resolvent analysis [16–18] has become widely used to model coherent structures in turbulent flows since it considers
the response of the linearised system to a forcing interpreted as the unknown non-linear term [19]. By singular
value decomposition (SVD) of the resolvent operator, optimal harmonic forcing modes and associated responses are
found. Resolvent analysis is used for the modelling of dominant coherent structures in turbulent flows [14, 19], data
assimilation [20–25], as well as flow control [26]. The comparison between resolvent analysis and SPOD has been
performed in turbulent channel flows at $Re_τ = 180$ and 550 in Abreu et al. [27], where good agreement has been
observed for elongated near-wall structures where the lift-up mechanism is active, associated with a dominance of the
first SPOD mode. The main limitation of the method lies in the fact that coherent structures are well predicted if
the non-linear term can be approximated as a Gaussian white noise, or if there is dominance of the first resolvent
amplification gain (singular value) [28]. These conditions are often not verified.

In the context of a triple decomposition, where the velocity field is split into a time-average, a coherent-structure
component and an incoherent turbulent field, an eddy viscosity can be introduced to the generalised Reynolds
stresses induced by the incoherent part [29]. For streaky structures in turbulent channel flows Cess’s eddy viscosity
model [30] has proven its prediction efficiency to some extent [24, 31–33], with a particular need in the logarithmic
layer. In Morra et al. [34], it has been shown that the cross spectral density (CSD) matrix of the non-linear forcing
in the DNS projects similarly onto resolvent forcing modes with eddy viscosity, thus explaining the improvement
by adding an eddy viscosity. In Amaral et al. [24], resolvent-based estimations have been performed in turbulent
channel flows at $Re_τ = 180, 550$ and 1000. It has been shown that in the buffer layer, both resolvent with and
without eddy viscosity lead to good estimations. However in the logarithmic layer, adding eddy viscosity becomes
necessary. Further improvement of estimation is possible if forcing statistics are used, without eddy viscosity, to
build an optimal estimator. This allows estimating flow fluctuations from wall measurements at various wall-normal
locations. Although using the CSD of the forcing terms is an interesting method to construct estimators, this is not a
viable approach to predict dominant coherent structures in turbulent flows, as it requires extensive use of flow data.

Resolvent analysis with Cess’s eddy viscosity will thus constitute our comparison model and will be referred to as
$ν_t$-resolvent analysis, in contrast with $ν$-resolvent analysis when no eddy viscosity is considered. This works well for
coherent structures, or waves, where strong production occurs. However, as argued by Symon et al. [35], since eddy
viscosity is mainly diffusive (up to eddy-diffusion gradients [36]), it breaks the energy conservation over the whole
spectrum. Then, it is not well adapted for waves receiving energy from other scales by backscattering. A detailed
study of the discrepancy of $ν$ and $ν_t$-resolvent analysis for a turbulent channel flow in terms of low-rank property,
projection onto SPOD modes and energy transfers can be found in Symon et al. [36]. Other attempts to improve
the modelling have been proposed. The embedding of covariance informations of the forcing has been for instance
proposed in [34, 37]. However, this strategy was considered for diagnostic purposes only and has not been considered
for predictions since a fine knowledge of the non-linear term is in that case required. An estimator has been proposed
by Gupta et al. [38] considering together eddy diffusion and a model of stochastic forcing. As an alternative in the
temporal domain, Zare et al. [39, 40] have devised a stochastic modelling based on control theory, which incorporates
a coloured-in-time noise.

In Tissot et al. [41], a modelling strategy based on stochastic transport, so-called stochastic linear modes (SLM),
has been proposed and will be considered in the present paper. It starts from a stochastic version of the Navier–Stokes
equations, originally introduced by Mémin [42], which is based on the stochastic transport of conserved quantities.
The formalism has been successfully employed to perform large eddy simulations [43], geophysical flow modelling [44–
49], near-wall flow modelling [50], data assimilation [51–53] and reduced-order modelling [54, 55]. An advantage of the
approach is the formulation of closure by defining statistics of a stochastic unresolved time-decorrelated (with respect
to the time scales of the contribution resolved by the model) velocity field. The associated random perturbation
ensues then from a stochastic transport operator. This stochastic transport involves in addition a stochastic diffusion,
and an effective drift velocity similar to the turbophoresis effect [56]. An exact energy balance is obtained between
the stochastic diffusion and the energy backscattering induced by the stochastic transport [44]. Linearising this model
and expressing it in the Fourier domain leads to what we refer to as *stochastic linear modes* (SLM).

In SLM, the non-linear term, interpreted as the wave-wave interactions, has been neglected, relying on the stochastic
transport of the solution by the incoherent turbulence to obtain a physically relevant model. In the present paper,
we come back to this strong assumption. The generation of streamwise vortices likely involves non-linear interactions
between large-scale coherent structures. As will be detailed further, a close analysis of SLM for these elongated
structures shows a poor prediction of the rolls, despite a good prediction of the streamwise velocity fluctuations. This
is consistent with the fact that coherent wave-wave interactions are neglected in SLM. In order to recover the right
roll properties, we propose in this paper to study the response of SLM to a “non-linear” forcing similarly to what is
done in resolvent analysis for modelling non-linear effects through an input-output formalism. We name this enhanced
solution forced stochastic linear modes (FSLM).

In addition to adding the aforementioned forcing, we propose some enhancements of the noise definition compared to Tissot et al. [41]. We propose an iterative procedure enforcing the noise to be incoherent with the solution. Moreover, the stochastic diffusion tensor is defined by root-mean-square velocity profiles, in order to ensure an approximated consistency between stochastic diffusion and noise expressed in the Fourier domain. We propose as well a decorrelation time definition based on an inertial range scaling. Finally, SLM/FSLM numerical computation is improved by the reformulation of the equations as an SVD problem.

With this more complete model which incorporates the effect of time-correlated turbulence on the coherent structures, we will explore the prediction abilities of stochastic modelling in the buffer and logarithmic layer of three turbulent channel flows at friction Reynolds number $Re_τ = 180$, $Re_τ = 550$ and $Re_τ = 1000$. In particular, we will explore the ability of FSLM to predict coherent structures in the logarithmic layer.

In section II, notations used along the paper are introduced. In section III, we present the stochastic model. In section IV we explore the ability of these models to predict buffer and logarithmic layer structures in turbulent channel flows. Some modelling recommendations are given in section V. Conclusions are provided in section VI. Presentation of the resolvent analysis, numerical details and complementary results are given in Supplementary Material in order to have a more complete view by varying Reynolds number and sweeping the wave-number space.

## II. NOTATIONS AND PRELIMINARIES

We consider three turbulent channel flows at the friction Reynolds numbers $Re_τ = 180$, $Re_τ = 550$ and $Re_τ = 1000$ with the Cartesian coordinates $\mathbf{x} = (x, y, z)$ of the streamwise, wall-normal and spanwise directions of the domain $\Omega$, respectively. The domain sizes ($L_x, L_y, L_z$) in outer units are respectively $(4\pi, 2, 2\pi)$, $(2\pi, 2, \pi)$ and $(2\pi, 2, \pi)$. Details and validations of the flow simulations can be found in [24], and additional details at $Re_τ = 550$ are present in [34]. The time-dependent ($t$) state variable $\mathbf{q}(x, y, z, t) = (\mathbf{u}, p)^T$ is composed of the velocity vector $\mathbf{u} = (u, v, w)^T$ and the pressure $p$. The velocity field is decomposed in its time-average and fluctuation $\mathbf{u} = \bar{\mathbf{u}} + \mathbf{u}'$ with $\bar{\mathbf{u}} = (U(y), 0, 0)^T$. By periodicity in the streamwise ($x$) and spanwise ($z$) directions, the space-time Fourier coefficient of the state variable with the sign convention $e^{i(\alpha x + \beta z - \omega t)}$ is noted $\hat{\mathbf{q}}_{\alpha, \beta, \omega}(y)$. Variables $\alpha, \beta, \omega$ refer respectively to streamwise wavenumber, spanwise wavenumber and angular frequency. In the wall-normal direction, we define a diagonal matrix $\mathbf{W}$ of quadrature coefficients. Finally, we note $\cdot^H$ the transpose-conjugate operation.

## III. STOCHASTIC LINEAR MODEL AND NON-LINEAR FORCING

### III.1. Stochastic linear modes

In Tissot et al. [41], a modelling strategy for coherent structures in turbulent flows has been proposed. The formalism relies on the stochastic transport of conserved quantities by a time-differentiable velocity component perturbed by the variation of a Brownian motion. Under these assumptions a stochastic version of the Navier-Stokes equations under location uncertainty [42] can be written. In this section, we recall how the stochastic model can be expressed in the frequency-wavenumber domain to predict coherent structures. More details can be found in Tissot et al. [41].

The displacement $\mathbf{X}(\mathbf{x}, t)$ of a particle is written in a differential form

$$d\mathbf{X}(\mathbf{x}, t) = \mathbf{u}(\mathbf{x}, t)dt + (\boldsymbol{\sigma}d\mathbf{B}_t)(\mathbf{x}), \quad (1)$$

where $\mathbf{u}$ is a time-differentiable velocity component, and $d\mathbf{B}_t$ is the increment of a Brownian motion. It can be remarked that equation (1) has to be understood as a time integral over an infinitesimal time increment $dt$. The operator $\boldsymbol{\sigma}$ is an integral operator which hides a spatial convolution in the domain $\Omega$ with a user-defined kernel $\tilde{\boldsymbol{\sigma}}$

$$(\boldsymbol{\sigma}d\mathbf{B}_t)^i(\mathbf{x}) = \int_{\Omega} \tilde{\boldsymbol{\sigma}}^{ij}(\mathbf{x}, \mathbf{x}', t)d\mathbf{B}_t^j(\mathbf{x}') d\mathbf{x}', \quad (2)$$

where the indices $i$ and $j$ refer to component indices. The vector $\mathbf{x}'$ is composed of the integration space coordinates.

Defined as such, $\boldsymbol{\sigma}d\mathbf{B}_t$ is the displacement field induced by a velocity component that is smooth in space, but decorrelated in time. This term aims at representing a time decorrelated (with respect to the time scale of the considered physical processes) turbulent velocity component. In the general framework, $\boldsymbol{\sigma}$ can be smoothly time-dependent, but for the application of the present paper in statistically stationary turbulent flows, we will assume it constant in time.
Associated with $\sigma$, we define the variance tensor $\mathbf{a}$ such that

$$ \mathbf{a}_{ij}(\mathbf{x})dt = \mathbb{E}((\sigma \mathrm{d}\mathbf{B}_t)^i(\mathbf{x}) (\sigma \mathrm{d}\mathbf{B}_t)^j(\mathbf{x})), \quad (3) $$

with $(\sigma \mathrm{d}\mathbf{B}_t)^i(\mathbf{x})$ the $i$th component of $\sigma \mathrm{d}\mathbf{B}_t$ at position $\mathbf{x}$ and $\mathbb{E}$ the expectation operator.

Using the Itō-Wentzell formula, conservation of mass and momentum subject to a stochastic transport leads to a stochastic version of the incompressible Navier-Stokes equations [42, 44], referred to as under location uncertainty:

$$ \begin{aligned} \mathrm{d}_t \mathbf{u} + (\mathbf{u}_d \cdot \nabla) \mathbf{u} \mathrm{d}t + (\sigma \mathrm{d}\mathbf{B}_t \cdot \nabla) \mathbf{u} = & -\nabla (p_t \mathrm{d}t + \mathrm{d}p_t) + \frac{1}{Re} \nabla \cdot (\nabla \mathbf{u}) \mathrm{d}t + \nabla \cdot \left(\frac{1}{2} \mathbf{a} \nabla \mathbf{u}\right) \mathrm{d}t + \frac{1}{Re} \nabla \cdot (\nabla \sigma \mathrm{d}\mathbf{B}_t) \\ \nabla \cdot \mathbf{u}_d = 0; \quad \nabla \cdot \sigma = 0, \end{aligned} \quad (4) $$

$$ \mathbf{u}_d = \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a}. $$

In system (4), $Re$ is the Reynolds number. Compared to the deterministic case, the transport of $\mathbf{u}$ by $\sigma \mathrm{d}\mathbf{B}_t$ is introduced. This term brings energy (backscatter) to the system, which is exactly compensated by the stochastic diffusion $\nabla \cdot (\frac{1}{2} \mathbf{a} \nabla \mathbf{u}) \mathrm{d}t$ [44]. The variable $\mathbf{u}_d$ is called drift velocity. It takes into account that, in average, particles tend to be transported from highly turbulent regions towards low-turbulence regions [see 54, and references therein]. Mass conservation leads to a divergence-free condition on $\mathbf{u}_d$, and on $\sigma$. Finally, a random pressure term $\mathrm{d}p_t$ corresponding to the small-scale velocity component is involved. This force balances the martingale part (proportional to $\mathrm{d}\mathbf{B}_t$) of the system.

In Tissot et al. [41], system (4) is linearised around a mean velocity profile $U(y)$ and written in the Fourier domain (see details of the derivation in the latter reference), leading to

$$ \begin{aligned} & -i\omega\hat{u}_{\alpha,\beta,\omega} + i\alpha U_d \hat{u}_{\alpha,\beta,\omega} + \hat{v}_{\alpha,\beta,\omega} \frac{\partial U}{\partial y} + i\alpha \hat{p}_{\alpha,\beta,\omega} + \tilde{D}(\hat{u}_{\alpha,\beta,\omega}) = -(\dot{\xi}_{\alpha,\beta,\omega})_y \frac{\partial U}{\partial y} + \frac{1}{Re} \Delta(\dot{\xi}_{\alpha,\beta,\omega})_x \\ & -i\omega \hat{v}_{\alpha,\beta,\omega} + i\alpha U_d \hat{v}_{\alpha,\beta,\omega} + \frac{\partial \hat{p}_{\alpha,\beta,\omega}}{\partial y} + \tilde{D}(\hat{v}_{\alpha,\beta,\omega}) = \frac{1}{Re} \Delta(\dot{\xi}_{\alpha,\beta,\omega})_y \\ & -i\omega \hat{w}_{\alpha,\beta,\omega} + i\alpha U_d \hat{w}_{\alpha,\beta,\omega} + i\beta \hat{p}_{\alpha,\beta,\omega} + \tilde{D}(\hat{w}_{\alpha,\beta,\omega}) = \frac{1}{Re} \Delta(\dot{\xi}_{\alpha,\beta,\omega})_z \\ & i\alpha \hat{u}_{\alpha,\beta,\omega} + \frac{\partial \hat{v}_{\alpha,\beta,\omega}}{\partial y} + i\beta \hat{w}_{\alpha,\beta,\omega} = 0 \quad ; \quad \frac{\partial \sigma_{xy}}{\partial y} = \frac{\partial \sigma_{yy}}{\partial y} = \frac{\partial \sigma_{zy}}{\partial y} = 0, \end{aligned} \quad (5) $$

with the modified diffusion operator

$$ \begin{aligned} \tilde{D}(\cdot) = & -\frac{1}{Re} \left(-\alpha^2 + \frac{\partial^2 \cdot}{\partial y^2} - \beta^2\right) \\ & -\frac{1}{2} \left(-\alpha^2 a_{xx} + i\alpha a_{xy} \frac{\partial \cdot}{\partial y} - \alpha\beta a_{xz} + i\alpha \frac{\partial a_{yx}}{\partial y} + \frac{\partial}{\partial y} \left(a_{yy} \frac{\partial \cdot}{\partial y}\right) + i\beta \frac{\partial a_{yz}}{\partial y} - \alpha\beta a_{zx} + i\beta a_{zy} \frac{\partial \cdot}{\partial y} - \beta^2 a_{zz}\right). \end{aligned} \quad (6) $$

The drift mean flow is $U_d(y) = U(y) - \frac{1}{2} \frac{\partial a_{xy}}{\partial y}$. The Fourier transform of $\sigma \mathrm{d}\mathbf{B}_t$ is noted $\mathrm{d}\xi_{\alpha,\beta,\omega}$, and the associated velocity Fourier component $\dot{\xi}_{\alpha,\beta,\omega} = \mathrm{d}\xi_{\alpha,\beta,\omega}/\mathrm{d}t$ is a standard centered Gaussian white noise convolved with the space-Fourier transform of $\sigma$. As the mean flow is parallel, the random transport term of equation (4) reduces to $-(\dot{\xi}_{\alpha,\beta,\omega})_y \frac{\partial U}{\partial y}$ in the right-hand-side of equation (5), which is the strain induced by extraction of energy of the mean flow by the turbulence. This term is central in the lift-up mechanism [7] and it is the main actor in the role of incoherent turbulence in the streaks of streamwise velocity $\mathbf{u}$. The choices of stochastic parameters $(\sigma, \mathbf{a})$ are detailed in section III.4.

The main added-value of SLM in wall-bounded flows is to model the impact of time-decorrelated turbulence on the lift-up mechanism and its associated momentum mixing by stochastic diffusion. Due to stochastic transport, equation (5) is a stochastic equation for $\hat{u}_{\alpha,\beta,\omega}$, Fourier transform of $\mathbf{u}'$. As a consequence, $\hat{u}_{\alpha,\beta,\omega}$ is a random variable, whose variability will allow us to extract purely coherent components through the estimation of the CSD matrix $\mathbb{E}(\hat{u}_{\alpha,\beta,\omega} \hat{u}_{\alpha,\beta,\omega}^*)$ and its eigenvectors. The leading eigenvector is called stochastic linear mode, or SLM. Our objective is to extract the dominant coherent component by SLM, which is compared to the leading spectral proper orthogonal decomposition (SPOD) mode [57]. In Tissot et al. [41], ensemble method is employed for the estimation and we present in section III.3 a reformulation of the problem as a singular value decomposition, to improve computational efficiency.
### III.2. Interactions between coherent structures

Equation (5) ensues from linearisation of equation (4). Coming back to the ground assumptions in our stochastic modelling, a triple decomposition is performed on the displacement

$$d\mathbf{X}(\mathbf{x}, t) = \bar{u}(\mathbf{x})dt + \mathbf{u}'(\mathbf{x}, t)dt + (\boldsymbol{\sigma}d\mathbf{B}_t)(\mathbf{x}). \quad (7)$$

The first term $\bar{u}(\mathbf{x})dt = (U(y) \ 0 \ 0)^T dt$ is the time-average displacement. The fluctuation is split in a time-differentiable component and an incoherent turbulent field, perceived as time-decorrelated compared to the time scale of the coherent structure, and modelled by a Brownian motion. Even if the time-average is non-ambiguous, the splitting of the fluctuation is less obvious in general and is often performed through a phase/ensemble average operator [29, 58]. The formulation (7) is a way to perform the triple decomposition in a unique manner through time differentiability of the variable (more precisely this decomposition is unique through the Bichteler-Dellacherie decomposition of stochastic processes [59]). Let us remark that the Brownian part $\boldsymbol{\sigma}d\mathbf{B}_t$ is modelled, while the time-differentiable part $\mathbf{u}'$ is solution of the system. Moreover, contrary to a splitting based on phase-averaging, $\mathbf{u}'$ contains coherent and incoherent contributions.

With this decomposition in mind, it can be seen that the stochastic diffusion can be interpreted as a generalised eddy diffusion (since a full tensor $\mathbf{a}$ is involved) induced by the noise. In this case, the diffusion comes directly from the time decorrelation assumption and stems from the Itō-Wentzel formula, where Itō quadratic variations can be viewed as providing local averaging coefficients. The diffusion does not come from a Boussinesq hypothesis. This diffusion term accounts for the effect of time-decorrelated component, and not for nonlinear interactions between time-differentiable components.

In system (5), the neglected term (written as a right-hand-side term in the momentum equation) is

$$\mathcal{F}\left(\left(\bar{\mathbf{u}}' \cdot \nabla\right)\bar{\mathbf{u}}' - \left(\mathbf{u}' \cdot \nabla\right)\mathbf{u}'\right), \quad (8)$$

where $\mathcal{F}(\cdot)$ stands for space-time Fourier transform. As in resolvent analysis (presented in Supplementary Material), this term is a convolution over all frequencies and wavenumbers, which renders an explicit expression difficult to obtain. A major difference compared to $\nu$-resolvent analysis is that it represents non-linear interactions between smooth-in-time structures carrying coherent wave contributions and does not include time-decorrelated turbulent fluctuations. In that sense, its interpretation is closer to the forcing term in $\nu_t$-resolvent analysis. We call the term (8) wave-wave interactions. The contribution of turbulent noise is already taken into account in the stochastic formulation.

We propose to treat the term (8) similarly to resolvent analysis, and to model it as a Gaussian white noise forcing term. The addition of a forcing term to equation (5) leads to

$$\left( \begin{array}{cccc} -i\omega + i\alpha U_d + \tilde{D}(\cdot) & \frac{\partial U}{\partial y} & 0 & i\alpha \\ 0 & -i\omega + i\alpha U_d + \tilde{D}(\cdot) & 0 & \frac{\partial \cdot}{\partial y} \\ 0 & 0 & -i\omega + i\alpha U_d + \tilde{D}(\cdot) & i\beta \\ i\alpha & \frac{\partial \cdot}{\partial y} & i\beta & 0 \end{array} \right) \left( \begin{array}{c} \hat{u}_{\alpha, \beta, \omega} \\ \hat{v}_{\alpha, \beta, \omega} \\ \hat{w}_{\alpha, \beta, \omega} \\ \hat{p}_{\alpha, \beta, \omega} \end{array} \right) = \left( \begin{array}{c} -(\dot{\xi}_{\alpha, \beta, \omega})_y \frac{\partial U}{\partial y} + \frac{1}{R_c} \Delta(\dot{\xi}_{\alpha, \beta, \omega})_x \\ \frac{1}{R_c} \Delta(\dot{\xi}_{\alpha, \beta, \omega})_y \\ \frac{1}{R_c} \Delta(\dot{\xi}_{\alpha, \beta, \omega})_z \\ 0 \end{array} \right) + b(y) \left( \begin{array}{c} \tilde{f}_x^{NL} \\ \tilde{f}_y^{NL} \\ \tilde{f}_z^{NL} \\ 0 \end{array} \right). \quad (9)$$

The linear operator in the left hand side of (9) can be written $\tilde{A}_{\alpha, \beta, \bar{q}} - i\omega\mathbf{E}$. The parameter $b(y)$ is an amplitude parameter of the non-linear forcings whose choice is based on the turbulent fluctuation level observed in the data. Its choice is described in section III.4. The vector $(\tilde{f}_x^{NL}, \tilde{f}_y^{NL}, \tilde{f}_z^{NL})^T$ carries independent standard centered Gaussian white noises. The above model will lead to forced stochastic linear modes, referred to as FSLM.

In the system (9), two stochastic right-hand side terms come from distinct physical mechanisms: the first term function of $\dot{\xi}_{\alpha, \beta, \omega}$ is related to stochastic transport by incoherent small scale turbulence, while the second forcing term accounts for the non-linear interactions between coherent structures.

### III.3. Numerical computation of the forced stochastic linear modes formalism

In Tissot et al. [41], an ensemble of solutions are computed to obtain an empirical CSD matrix. This procedure turns out to be more expensive than a singular value decomposition (SVD) for small size problems (one-dimensional in the
y direction). We can note that for large scale problems, advanced ensemble-based techniques [60, 61] or time-domain formulations [62] can be employed. We propose here to write FSLM as an SVD problem. Starting from system (9) and similarly as in resolvent analysis (presented in supplementary material), we define

$$ \tilde{\mathbf{L}}_{\alpha, \beta, \omega} = \mathbf{H} (\tilde{\mathbf{A}}_{\alpha, \beta, \bar{q}} - i\omega \mathbf{E})^{-1} \tilde{\mathbf{B}}, \quad (10) $$

with

$$ \mathbf{H} = \begin{pmatrix} \mathbb{I} & 0 & 0 & 0 \\ 0 & \mathbb{I} & 0 & 0 \\ 0 & 0 & \mathbb{I} & 0 \end{pmatrix}, \quad \tilde{\mathbf{B}} = \begin{pmatrix} -(\mathbf{\Phi}^{\sigma})_{y} \mathbf{D}^{\sigma} \frac{\partial \mathbf{U}}{\partial y} + \frac{1}{Re} \Delta(\mathbf{\Phi}^{\sigma})_{x} \mathbf{D}^{\sigma} & b(y)\mathbb{I} & 0 & 0 \\ \frac{1}{Re} \Delta(\mathbf{\Phi}^{\sigma})_{y} \mathbf{D}^{\sigma} & 0 & b(y)\mathbb{I} & 0 \\ \frac{1}{Re} \Delta(\mathbf{\Phi}^{\sigma})_{z} \mathbf{D}^{\sigma} & 0 & 0 & b(y)\mathbb{I} \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad (11) $$

with $\mathbb{I}$ the identity matrix. The output operator $\mathbf{H}$ specifies that the modes will be optimal in terms of kinetic energy. The input operator $\tilde{\mathbf{B}}$ maps the vector of random variables

$$ \tilde{\mathbf{f}} = (\eta_1, \dots, \eta_{N_{\sigma}}, \tilde{f}_{x}^{NL}, \tilde{f}_{y}^{NL}, \tilde{f}_{z}^{NL})^T $$

to the forcing space of the linearised system. The matrix $(\mathbf{\Phi}^{\sigma})_j$, for $j = \{x, y, z\}$, gathers in columns $(\mathbf{\Phi}_k^{\sigma})_j$ for $k \in [1, N_{\sigma}]$, an expansion basis of the noise $d\xi_{\alpha, \beta, \omega} = \sum_{k=1}^{N_{\sigma}} c_k \mathbf{\Phi}_k^{\sigma} \eta_k$ which will be specified in section III.4, with $\eta_k$ standard centered Gaussian white noises. The diagonal matrix $\mathbf{D}^{\sigma}$ contains the associated amplitude coefficients $c_k$.

We perform the singular value decomposition

$$ \mathbf{W}^{\frac{1}{2}} \mathbf{H} (\tilde{\mathbf{A}}_{\alpha, \beta, \bar{q}} - i\omega \mathbf{E})^{-1} \tilde{\mathbf{B}} \mathbf{W}_{f}^{-\frac{1}{2}} = \mathbf{U}_{\text{FSLM}} \mathbf{\Sigma}_{\text{FSLM}} \mathbf{V}_{\text{FSLM}}^*, \quad (12) $$

with

$$ \mathbf{W}_f = \begin{pmatrix} \mathbb{I} & 0 \\ 0 & \mathbf{W} \end{pmatrix}. \quad (13) $$

As a final step, FSLM are defined by $\mathbf{\Phi}_i^{\text{FSLM}} = \mathbf{W}^{-\frac{1}{2}} \mathbf{V}_{\text{FSLM},i}$, where the first mode is the predicted coherent structure and the higher order modes are used to define the noise at the next iteration (see section III.4). Moreover, as in resolvent analysis, an estimation of the CSD matrix can be obtained by $\mathbf{S} = \tilde{\mathbf{L}}_{\alpha, \beta, \omega} \tilde{\mathbf{L}}_{\alpha, \beta, \omega}^*$. 210 211 212

### III.4. Choice of parameters in FSLM

We recall that we focus on coherent structures perturbed by turbulent flows in the buffer and logarithmic layers at scales where production exceeds dissipation, for which, therefore, a forward energy cascade is expected [35]. In the logarithmic layer, we focus at energetic scales, and as highlighted in Jiménez [63], dissipation takes place at a smaller scale. We expect in this region an energy cascade draining energy from large to small scales through an inter-scale energy flux. For large energetic scales in the logarithmic layer, we expect as well a larger influence of incoherent turbulence onto the wave compared to the buffer layer; we aim at modelling such influence by FSLM.

The two-point statistics of the noise, carried by $\sigma$, have to represent time decorrelated turbulent velocity field fluctuations. Its definition is an open question and relies on an a priori knowledge of the fluctuating velocity field. Our strategy is to use few parameters, preferably with quantities available in standard simulation data or well documented in the literature. Moreover, we need to respect the ground hypothesis that the noise is decorrelated from resolved coherent field (at the large-scale characteristic time scale).

We propose to set the variance tensor, $\mathbf{a}$, defined in equation (3), from root-mean-square (RMS) velocity profiles, and variances of velocity fluctuations, which are quantities often available in databases accompanying the mean flow profile:

$$ \mathbf{a}(y) = \tau_s \begin{pmatrix} \langle u'^2(y) \rangle & \langle u'(y)v'(y) \rangle & 0 \\ \langle u'(y)v'(y) \rangle & \langle v'^2(y) \rangle & 0 \\ 0 & 0 & \langle w'^2(y) \rangle \end{pmatrix}, \quad (14) $$

where $\langle \cdot \rangle$ denotes average in time and in the homogeneous directions. The underlying hypothesis to use RMS profiles is that the contribution of the single coherent wave we are trying to predict is small compared to the whole time-domain solution. Thus, the RMS, which contains all contributions of the turbulent velocity field is a fair estimate of
the turbulence which impacts the wave. The decorrelation time $\tau_s$, necessary for dimensional consistency, represents the time scale necessary for the Brownian motion to perform mixing by stochastic diffusion. This parameter is crucial for obtaining relevant results since it controls the level of diffusion. The time scale $\tau_s$ should represent, at a given wavelength, the time scale necessary for the turbulence to affect the wave by a transport mechanism. For this, we rely on an inertial scaling $\tau_s = \tau_0(l/l_0)^{\frac{3}{2}}$ proposed in [64], assuming that the wavelength lies within the inertial range of an energy cascade under Kolmogorov hypotheses, and that this turbulent field is incoherent with the wave solution. The time $\tau_0 = l_0/U_0$ is the outer time scale, $l_0 = 2$ is the channel height, $U_0$ is the velocity averaged over the wall-normal direction; the scale of the wave is $l = 2\pi/\sqrt{k_x^2 + k_z^2}$ with $k_x = 2\pi/\lambda_x$ and $k_z = 2\pi/\lambda_z$ [10]. This scaling is valid for scales such that $l < l_0$. We do not expect our scaling to be valid for $l$ larger than the channel height leading to structures living in the outer region. It can be noticed that the structure of the model allows a scale-dependent stochastic diffusion through the decorrelation time $\tau_s$, which is set here by a physical scaling. In Gupta et al. [38], a similar scale dependence of the eddy diffusion has been observed to be necessary to produce accurate results.

The noise $d\xi_{\alpha,\beta,\omega}$ is the space-time Fourier transform of $\sigma dB_t$. It is white in time, and its covariance is the CSD of $\sigma dB_t$. It should match the Fourier transform of the cross correlation tensor of $\sigma dB_t$, whose diagonal is the variance tensor $\mathbf{a}$. Indeed, the cross spectral density of $\sigma dB_t$ is

$$ \mathbb{E} \left( d\xi_{\alpha,\beta,\omega}^i \left( d\xi_{\alpha,\beta,\omega}^j \right)^H \right) = \mathcal{F} \left( \mathbb{E} \left( \left( \sigma dB_t \right)_x^i \left( \sigma dB_{t'} \right)_{x'}^j \right) \right), \quad (15) $$

and we recall the link with the variance tensor $\mathbf{a}$ in equation (3). As a consequence, defining the noise based on an expansion onto SPOD modes $\Phi_k^{\text{SPOD}}$ associated with the eigenvalues $\lambda_k^{\text{SPOD}}$

$$ d\xi_{\alpha,\beta,\omega} = \sum_{k=1}^{\infty} \sqrt{\lambda_k^{\text{SPOD}}} \Phi_k^{\text{SPOD}} \eta_k, \quad (16) $$

with $\eta_k$ standard centered Gaussian white noise variables, is consistent with the definition of $\mathbf{a}$ in equation (14), since SPOD modes are eigenfunctions of the CSD matrix of the turbulent fluctuation. This definition has two main issues. The first is the requirement of a fine description of the turbulent velocity field, since SPOD modes are required. Instead, we aim at constructing a model with a reduced quantity of data. The second issue is that in definition (16), the FSLM coherent structure would be necessarily correlated with a part of the noise, since it would be spanned by the full SPOD basis. Ultimately, we would like to subtract it from the noise basis, but only at the considered frequency-wavenumber couple to not invalidate our definition (14).

To address these issues, we propose to relax the strong consistency between equation (16) and the diffusion tensor $\mathbf{a}$ equation (14). Instead, we ensure that the noise is decorrelated from the wave by modifying the definition of the noise. First, we express it as an expansion onto an orthonormal basis

$$ d\xi_{\alpha,\beta,\omega} = \sum_{k=1}^{N_{\sigma}} c_k \Phi_k^{\sigma} \eta_k. \quad (17) $$

We then propose a first guess by defining $\Phi_k^\sigma = \Phi_{k+1}^{\nu\text{-resolvent}}$ and $c_k = \sqrt{\lambda_1^{\text{SPOD}} s_{k+1}^{\nu\text{-resolvent}}/s_1^{\nu\text{-resolvent}}}$, with $k \in [1, \dots, N_\sigma]$, ($s_k^{\nu\text{-resolvent}}, \Phi_k^{\nu\text{-resolvent}}$) the $k^{\text{th}}$ singular value and optimal response mode of $\nu$-resolvent analysis and $\lambda_1^{\text{SPOD}}$ the first SPOD eigenvalue. This guess rescales the noise spanned by $\nu$-resolvent suboptimal modes in such a way that the energy of the first mode matches the first SPOD mode. Doing this, we define an orthonormal family of vectors orthogonal to the dominant resolvent mode. The amplitude rescaling aims at obtaining a rough approximate consistency between $\sigma$ and the definition of $\mathbf{a}$ by the RMS profiles. The use of resolvent modes as a first guess frees the modelling from the data. Only the first SPOD eigenvalue is required, but this single parameter can be replaced by a free parameter fixed by some physical knowledge to obtain a full model-based procedure, in which no data is required.

In a second step, we correct the definition of (17) in order to ensure that the noise is decorrelated from the first FSLM. Once expressed in the Fourier domain, the time decorrelation becomes a decorrelation between ensemble realisations of the Fourier component. As explained in Towne et al. [57] for the SPOD modes, since the CSD matrix has been diagonalised, the contribution of separate eigenfunctions are decorrelated. We then construct a noise spanned by the eigenfunctions of the CSD matrix excluding the first FSLM. For that, we choose $\Phi_k^\sigma = \Phi_{k+1}^{\text{FSLM}}$ and $c_k = \sqrt{\lambda_{k+1}^{\text{FSLM}}}$, with $k \in [1, \dots, N_\sigma]$, where ($\lambda_k^{\text{FSLM}}, \Phi_k^{\text{FSLM}}$) are eigen-elements of CSD matrix $\mathbf{S}$ of FSLM solutions. The procedure is cyclic since high-order modes of FSLM are mandatory to predict the leading FSLM mode, but it is possible to compute it iteratively through a fixed point procedure initialised with the first guess, as summarised in algorithm 1. In practice, calculations converge quickly in few (less than 10) iterations with a relative tolerance on the Frobenius norm $\| \cdot \|_F$ of the CSD equal to $\epsilon = 10^{-3}$. An example of convergence is shown in the Supplementary Material. It
can be noticed in particular that it converges toward a solution where the energy decay of the spectrum is similar to SPOD (see section IV.4), which does not invalidate the consistency between noise and stochastic diffusion. Moreover, the decay in the singular values is moderately fast suggesting an incoherent (by construction) contribution that is sufficiently small to be considered as a noise, but too large to be neglected.

<caption>**Algorithm 1:** Iterative procedure for FSLM</caption>

Compute $\nu$-resolvent: $\mapsto (S^{\nu\text{-resolvent}}, \Phi^{\nu\text{-resolvent}})$
$\Phi_k^{\sigma,(1)} \leftarrow \Phi_{k+1}^{\nu\text{-resolvent}}$, $c_k^{(1)} = s_{k+1}^{\nu\text{-resolvent}} \frac{\sqrt{\lambda_1^{\text{SPOD}}}}{s_1^{\nu\text{-resolvent}}}$, $S^{(1)} \leftarrow \sum_{k=0}^{N_\sigma} (c_k^{(1)})^2 \Phi_k^{\sigma,(1)} \Phi_k^{\sigma,(1)*}$;
$n \leftarrow 1$;
**while not converged do**
&nbsp;&nbsp;&nbsp;&nbsp;Compute FSLM by SVD procedure (sec. III.3): $\mapsto (\tilde{L}_{\alpha,\beta,\omega}, \lambda^{\text{FSLM}}, \Phi^{\text{FSLM}})$;
&nbsp;&nbsp;&nbsp;&nbsp;$S^{(n+1)} \leftarrow \tilde{L}_{\alpha,\beta,\omega} \tilde{L}_{\alpha,\beta,\omega}^*$;
&nbsp;&nbsp;&nbsp;&nbsp;$\Phi_k^{\sigma,(n+1)} \leftarrow \Phi_{k+1}^{\text{FSLM}}$, $c_k^{(n+1)} \leftarrow \sqrt{\lambda_{k+1}^{\text{FSLM}}}$;
&nbsp;&nbsp;&nbsp;&nbsp;if $\|S^{(n+1)} - S^{(n)}\|_F / \|S^{(1)}\|_F < \epsilon$ then
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;converged $\leftarrow True$;
&nbsp;&nbsp;&nbsp;&nbsp;end
end

To summarise, the proposed procedure uses the RMS profiles to define the diffusion tensor **a**, which is the one defined in the time-domain in equation (3). The noise $d\xi_{\alpha,\beta,\omega}$, space-time Fourier transform of $\sigma dB_t$, is expanded on an orthonormal basis, which is estimated by an iterative procedure ensuring decorrelation between the noise and the solution. An initial guess is defined by resolvent modes rescaled using the first SPOD eigenvalue in order to obtain consistency with the definition of **a** with a minimum of data.

Finally, in FSLM, the non-linear forcing amplitude has to be given. In order to obtain a physically relevant order of magnitude, the non-linear forcing amplitude $b(y)$ is chosen as $(\sqrt{\lambda_1^{\text{SPOD}}}/s_1^{\nu\text{-resolvent}}) (\text{TKE}(y)/\max(\text{TKE}))$, with $\text{TKE}(y) = \langle u'^2 \rangle + \langle v'^2 \rangle + \langle w'^2 \rangle$ the turbulent kinetic energy. This scaling allows us to define a profile of non-linear forcing in the wall-normal direction consistent with the turbulent activity, and such that the response of the deterministic linearised system (without eddy viscosity) to this forcing has an amplitude comparable with SPOD.

## IV. APPLICATION TO TURBULENT CHANNEL FLOW

### IV.1. Numerical simulation

Databases of direct numerical simulation of turbulent channels were obtained with the pseudospectral code Channelflow 2.0 [65]. Periodic boundary conditions are enforced in the streamwise ($x$) and spanwise ($z$) directions and Chebyshev polynomials are used in the wall-normal direction ($y$). Parameters are given in table I and additional numerical details, including validation results, can be found in Amaral et al. [24]. Mean flow profiles for the three Reynolds numbers and root-mean-square (RMS) profiles at $Re_\tau = 1000$ are presented in figure 1. Results are presented using non-dimensional quantities using viscous (wall) scaling, denoted with a + superscript.

SPOD has been computed as a reference to which $\nu_t$-resolvent analysis and FSLM will be compared. They represent the most energetic structure for a given frequency-wavenumber combination, which can be meaningfully compared to most amplified responses of resolvent and FSLM. Complementary numerical details are given in Supplementary Material.

Resolvent modes are known [19, 60] to show large responses around the critical layer $y_c^+$, i.e. where the phase speed $c^+ = \lambda_x^+ / \lambda_t^+$ matches the mean flow $U^+(y_c^+)$. SPOD modes follow the same trend, which is consistent with the fact that these modes are equivalent if the non-linear term behaves as a Gaussian white noise [66]. Two waves have been selected: one denoted W1 with $(\lambda_x^+, \lambda_z^+, \lambda_t^+) \approx (1000, 100, 100)$, typical of the streaks structures in the buffer layer

<table><caption>TABLE I. Numerical parameters for the simulations.</caption><thead><tr><th>Re<sub>τ</sub></th><th>Re</th><th>N<sub>x</sub></th><th>N<sub>y</sub></th><th>N<sub>z</sub></th><th>Δx<sup>+</sup></th><th>Δy<sub>min</sub><sup>+</sup></th><th>Δy<sub>max</sub><sup>+</sup></th><th>Δz<sup>+</sup></th><th>Δt<sup>+</sup></th></tr></thead><tbody><tr><td>180 (179)</td><td>2800</td><td>192</td><td>129</td><td>192</td><td>11.7</td><td>5.4 · 10<sup>-2</sup></td><td>4.4</td><td>5.9</td><td>5.7</td></tr><tr><td>550 (543)</td><td>10000</td><td>384</td><td>257</td><td>384</td><td>8.9</td><td>4.1 · 10<sup>-2</sup></td><td>6.7</td><td>4.4</td><td>3.0</td></tr><tr><td>1000 (996)</td><td>20000</td><td>484</td><td>385</td><td>484</td><td>12.9</td><td>3.3 · 10<sup>-2</sup></td><td>8.2</td><td>6.5</td><td>2.5</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Mean and root-mean-square profiles. Grey areas indicate the spatial supports of W1 and W3.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Effect of the drift velocity.</figcaption></figure>

chosen at the lowest Reynolds number $Re_τ = 180$; and one denoted W3 with $(\lambda_x^+, \lambda_z^+, \lambda_t^+) \approx (2000, 500, 100)$, evolving within the logarithmic layer at the highest Reynolds number $Re_τ = 1000$. In Supplementary Material, numerical details are presented. Moreover, the robustness of the method is shown by varying Reynolds number for W1, and two other waves (W2 and W4) are presented in order to vary the wall-distance of the wave spatial support. As in Tissot et al. [41] we consider modes that are odd in $u$ and $w$ (and thus even in $v$) around the channel centerline. This has been performed by enforcing symmetry in the operators $L_{\alpha,\beta\omega}$, $H$, $\bar{H}$.

The spatial supports of the waves W1 and W3 are reported in figures 1. In figure 2(a), the drift velocity (in wall units) associated with W1 to W4 are displayed at $Re_τ = 1000$. It shows that the corrective drift $-\frac{1}{2}\partial a_{xy}/\partial y$ plays essentially a role in the buffer region. Since the variance tensor $a$ is defined based on RMS profiles in equation (14), the drift velocity accounts for the effective transport induced by the wall-normal variations of $\langle u'(y)v'(y) \rangle$. In the buffer region, the magnitude of $\langle u'v' \rangle$ (with negative values) increases with the wall distance, and this inhomogeneity tends to induce a positive streamwise transport velocity, which corresponds to the corrective drift term $-\frac{1}{2}\partial a_{xy}/\partial y$. This effective transport by the turbulent fluctuations is taken into account in the proposed linearised model, and does not appear explicitly in an eddy-viscosity model. This relevance of the drift velocity in the buffer region is in line with the observations and modelling of Pinier et al. [50]. Additionally, we can see that with our definition, the effect of drift velocity is more pronounced for waves evolving at higher wall-normal distance due to larger decorrelation times $τ_s$. As a matter of fact, for such long waves there is a more substantial contribution of the stochastic transport, which occurs with a longer decorrelation time.
<figure><img src="image_4.png" /><figcaption>FIG. 3. Reconstructions of W1 at $Re_τ = 180$. Colors are streamwise velocities, arrows are in-plane velocity fields.</figcaption></figure>

IV.2. Buffer layer

In the buffer layer, we present the results at $Re_τ = 180$, and complementary results at $Re_τ = 550$ and $Re_τ = 1000$ are given in the Supplementary Material. Figure 3(a) shows a velocity field cross section of the SPOD for W1 ($λ_x^+ = 1124$, $λ_z^+ = 102$, $λ_t^+ = 100$). It shows a typical streaky structure of $u$ with streamwise vortices (rolls), which highlights the lift-up mechanism: in regions of high streamwise velocity high-speed streak are emerging. They are associated with negative $v$ components, which transport fluid with high streamwise velocity to a region with lower mean flow; the opposite happens for low-speed streaks, which are associated with positive $v$ components. Predictions by $v$-resolvent and $v_t$-resolvent analysis are shown in figures 3(b) and 3(c) respectively. We can see a relevant prediction, with an improvement when eddy viscosity is added. This is consistent with Morra et al. [32]. Figure 3(d) shows the solution of the proposed stochastic model, but omitting the non-linear forcing (SLM). It can be seen that the streaks are well predicted, but the rolls are absent. However, taking into account wave-wave interactions by a non-linear forcing (FSLM) enables us to recover the rolls, and to obtain accurate predictions. This can be explained by the fact that stochastic transport models the effect of the incoherent part of the velocity field, thus leading to good predictions of the $u$ profiles. Since near-wall streamwise vortices are thought to arise from a non-linear interaction of coherent structures [5], a non-linear forcing is mandatory to predict them. Resolvent analysis with eddy viscosity leads to good predictions as it takes into account this non-linear forcing and incorporates eddy diffusion. These predictions are significantly enhanced by the stochastic model since it explicitly modifies lift-up by incoherent turbulent motions through three complementary terms: the transport by the noise, a diffusion tensor with non-zero off-diagonal terms and a drift velocity active in the buffer region [41].

Profiles of power spectral density (PSD) of the three velocity components are shown in figure 4. They confirm that streamwise velocity ($u$) profiles are similarly captured by $v_t$-resolvent and by SLM. The agreement of FSLM with SPOD data is significantly improved. Moreover, the streamwise vortices signing on the ($v, w$) profiles are not captured by SLM but strongly intensified in FSLM. The wall-normal velocity ($v$) is especially affected by the stochastic transport leading to the best agreement with SPOD.

To demonstrate the robustness of the procedure, buffer layer modes at other Reynolds numbers are shown in Supplementary Material. Despite a slight overall deterioration of agreement between all models and SPOD when the Reynolds number increases, the trend is maintained and FSLM shows systematically a better agreement.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 4. PSD velocity profiles of W1 at $Re_τ = 180$.</figcaption></figure>

347

IV.3. Logarithmic layer

We now select a wave named W3 evolving within the logarithmic layer by setting the phase speed $c^+ = 18.1$ associated with a critical level $y_c^+ = 180$. The wave-number ($\lambda_x^+ = 2087$, $\lambda_z^+ = 522$) has been chosen in the energy peak [taken from ref. 24] of the premultiplied powerspectra of the streamwise velocity. It can be extracted by SPOD, as shown in figure 5(a). The $\nu$-resolvent analysis (figure 5(b)) extracts typical critical layer modes, with a narrow spatial support located at the critical layer, i.e. at the wall-normal position $y_c^+$ where the phase speed $c^+$ matches the mean velocity. Incorporating eddy viscosity (figure 5(c)) leads to wider spatial support, more similar to SPOD modes, which is again consistent with previous studies [32]. However, there is room for improvements, since SPOD modes show a structure that peaks further from the wall than what is predicted by $\nu_t$-resolvent. FSLM in figure 5(d) improves significantly this prediction with the streamwise velocity structure further from the wall and a more accurate shape of the rolls.

348

349

350

351

352

353

354

355

356

357

Figure 6 showing the PSD profiles highlights quantitatively this improvement. Concerning the $u$ component, figure 6(a) shows that differently from the other models, the spatial support is very well captured. As for the wall-normal $v$ velocity, the shape of the profile is better predicted but with a high relative amplitude. The spanwise $w$ velocity is better captured as well.

358

359

360

361

The lower accuracy of $\nu_t$-resolvent predictions can be understood by the fact that the effect of the incoherent turbulent field on the wave is modelled only by a diffusive mechanism. On the contrary, FSLM incorporates through stochastic transport some driving mechanisms induced by the incoherent motions existing at the same scale. The success of FSLM suggests that in the logarithmic region, where the turbulence is developed, taking into account the stochastic nature of the log-layer structures is central to perform accurate predictions.

362

363

364

365

366

367

In addition, profiles of SLM, i.e. neglecting the non-linear forcing, are shown to produce poor predictions. This suggests again that coherent non-linear wave-wave interactions are crucial for self-sustaining process for log-layer. It corroborates hypotheses in Flores and Jiménez [8] and Cossu and Hwang [12] that a coherent large scale self-sustaining

368

369
<figure><img src="image_8.png" /><figcaption>FIG. 5. Reconstructions of W3 at Re$_{τ}$ = 1000.</figcaption></figure>

372 process is in action for large log-layer structures.
Frequency-wavenumber space has been swept, and colinearity metrics

$$
\theta_{\alpha, \beta, \omega}^{\text{model}} = \frac{\left| \left( \mathbf{\Phi}_{1, \alpha, \beta, \omega}^{\text{model}}, \mathbf{\Phi}_{1, \alpha, \beta, \omega}^{\text{SPOD}} \right) \right|}{\left| \mathbf{\Phi}_{1, \alpha, \beta, \omega}^{\text{model}} \right| \left| \mathbf{\Phi}_{1, \alpha, \beta, \omega}^{\text{SPOD}} \right|} \quad (18)
$$

373 have been computed (metric used for instance in [67]). This metric is a normalised inner-product between the dominant
374 SPOD mode $\mathbf{\Phi}_{1,\alpha,\beta,\omega}^{\text{SPOD}}$ and the mode issued from a given model $\mathbf{\Phi}_{1,\alpha,\beta,\omega}^{\text{model}}$. A value of 1 means exact collinearity between
375 modes, while 0 happens when the modes are orthogonal. Then, we compute the metric $\gamma_{\alpha,\beta,\omega} = \log (\theta_{\alpha,\beta,\omega}^{\text{FSLM}} / \theta_{\alpha,\beta,\omega}^{\nu t\text{-resolvent}})$
376 which represents the improvement ($\gamma > 0$) or deterioration ($\gamma < 0$) of colinearity with SPOD compared to the $\nu t$-
377 resolvent model. Figure 7 shows the value of $\gamma_{\alpha,\beta,\omega}$ at four critical layer positions as a function of streamwise and
378 spanwise wavenumbers. We can see that in the buffer and logarithmic layer, a wide range of streamwise elongated
379 structures are improved with FSLM compared to $\nu t$-resolvent analysis. The improvement is more pronounced further
380 from the wall since agreement is more difficult to obtain. Complementary maps of $\theta_{\alpha,\beta,\omega}^{\text{model}}$ are given in Supplementary
381 Material. Isocontour of the pre-multiplied first SPOD eigenvalue $\alpha\beta\lambda_1^{\text{SPOD}}$ are superimposed, and show that deteri-
382 oration happens at scales where less energy is present. Finally, in figure 7(d) almost in the outer-region, we can see
383 that FSLM provide slightly worse performances than $\nu t$-resolvent for large $\lambda_x$ and $\lambda_z$. We explain this discrepancy by
<figure><img src="image_9.png" /><figcaption>(a) Streamwise velocity $|\hat{u}|^2$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Wall-normal velocity $|\hat{v}|^2$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Spanwise velocity $|\hat{w}|^2$.</figcaption></figure>

<caption>FIG. 6. PSD velocity profiles of W3 at $Re_τ = 1000$.</caption>

the choice of decorrelation time $\tau_s$ which is designed based on inertial-range scalings (see section III.4). These maps prove a wide range of validity of the proposed model.

IV.4. Eigenspectrum

In figure 8, we compare the eigenvalues $\lambda^{\text{FSLM}}$ of the FSLM model, i.e. the eigenvalues of $S^{(n)}$ once the iterative procedure is converged, with the SPOD eigenvalues. It is performed for the two examples W1 at $Re_τ = 550$ and W3 at $Re_τ = 1000$. As a reference, we show the spectrum of $\nu$-resolvent amplification energy gain $(\sigma^{\nu\text{-resolvent}})^2$, which constitutes the initial guess of the iterative procedure. It can be seen that even if the initial guess has a spectrum with a decay that is too fast, the iterative procedure succeeds to produce a rate of energy decay similar to SPOD. We recall that these values have not been informed in the model, and it is only a consequence of the constraint that the noise is decorrelated from the FSLM coherent structure. We recall that a noise based on an expansion onto SPOD ensures consistency between the noise expressed in the Fourier domain and the definition given to the stochastic diffusion using RMS profiles. The fact that we recover a spectrum similar to SPOD suggests that the incoherent field spanned by suboptimal is relevant, despite the relaxation of the strong consistency between SPOD and RMS profile. Moreover, the moderate decay in the spectrum indicates that, after convergence, the incoherent part cannot be neglected.

As an indication, we have also shown in figure 8 the energy gain spectrum of $\nu_t$-resolvent analysis. Consistently with the observations in [36], in the energetic scales of the buffer layer the energy gain of suboptimal modes of $\nu_t$-resolvent analysis decays as fast as in $\nu$-resolvent analysis, with a strongly low-rank behaviour. However in the log-layer, adding eddy-viscosity reduces this low-rank behaviour, and we can see in figure (8) that it rejoins the decay rate of SPOD. FSLM has a decay rate of the spectrum close to SPOD in both cases.

As a caveat, we recall that the comparison is made with a specific eddy viscosity model, which has been adjusted in Del Álamo and Jiménez [68] to match the numerical mean flow at $Re_τ = 2000$. It may be highlighted that this eddy viscosity is not necessarily optimal for all Reynolds numbers. Moreover, the eddy viscosity designed to match
<figure><img src="image_12.png" /><figcaption>FIG. 7. Metric $\gamma_{\alpha,\beta,\omega} = \log(\theta_{FSLM}/\theta_{\nu_t\text{-res}})$ of improvement ($> 0$) or deterioration ($< 0$) of collinearity between FSLM and SPOD compared to $\nu_t$-resolvent analysis as a function of $\lambda_x$, $\lambda_z$ for various critical layer positions. Isocontours are the pre-multiplied value of the first SPOD $\alpha\beta\lambda_1^{\text{SPOD}}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 8. Comparison of the eigenvalues. SPOD eigenvalues $\lambda^{\text{SPOD}}$, $\nu$-resolvent amplification energy gain $(\sigma^{\nu\text{-resolvent}})^2$ and FSLM spectrum $\lambda^{\text{FSLM}}$ (eigenvalues of $S^{(n)}$ once converged). The spectra have been normalised by the first SPOD eigenvalue.</figcaption></figure>

the mean flow is also not necessarily optimal to predict a coherent structure at a given wavenumber/frequency couple. There are, indeed, some indications in the literature that the eddy viscosity should be dependent on the wavenumber [36, 38]. On the other hand, the FSLM incorporates as well a given amount of data, but no optimisation procedure to match data. A fully consistent comparison is not obvious, and we keep in mind that some amount of improvement could be obtained as well in $\nu_t$-resolvent analysis by optimising an eddy viscosity parameter. The goal is here more to highlight the benefit of stochastic modelling under location uncertainty than putting models in competition.

V. MODELLING RECOMMENDATIONS

The proposed formalism brings some modelling tools which address some limitations of the resolvent analysis framework. The non-linear term, assumed to act as an additive Gaussian white noise in resolvent analysis, is constrained in SLM to be issued from a transport mechanism. This can be interpreted in the resolvent framework as playing the role of colouring the additive forcing term. Actually it is more than this, since it takes explicitly into account inho-
mogeneity and anisotropy of the turbulent fluctuations. Besides, as explored for instance in Gupta et al. [38], there is often – on the basis for instance of the fluctuation dissipation theorem – the requirement to model a random term which brings energy to the system and a diffusion term which dissipates it. The stochastic framework under location uncertainty provides by construction an energy balance between these two mechanisms. For all these reasons, despite its more complex structure, stochastic modelling under location uncertainty may be relevant in flow configurations where resolvent analysis requires finer modelling.

In channel flows, $\nu$-resolvent is accurate only close to the wall for elongated structures (see supplementary material for an overview of the metric $\theta_{\alpha,\beta\gamma}$). Adding eddy viscosity becomes necessary at higher wall-normal locations. In general flow configurations, we do not have systematically a model of eddy viscosity. In such a case a Reynolds Averaged Navier-Stokes (RANS) calculation, which requires a closure as well, allows to obtain eddy diffusion, as performed for instance in Pickering et al. [69]. In the proposed formalism, the stochastic diffusion is deduced from some knowledge of the small scale statistics, expressed conveniently in terms of velocity fluctuations.

In the present paper, we add the contribution of time-coherent non-linear interactions. In flows where the main effect lies in small-scale turbulence, the forcing term may be unnecessary. Dedicated studies should be performed to answer this question. However, the non-linear forcing term appears to be important when non-linearities arise from interactions between time-coherent structures, such as in wall-bounded flows.

SLM and FSLM are characterised by a modelling degree of freedom through the choice of the parameters ($\bar{u}$, $\sigma$, $\mathbf{a}$, $\tau_s$). We note in particular the role of $\tau_s$, which may be sensitive since it controls the amplitude of the stochastic diffusion. These parameters can be determined based on data, or some physical scalings. The proposed strategy in section III.4 is not universal, but aims at minimising the data requirement. Another possible strategy, which is currently being pursued by our group, is to obtain these parameters without data, but only based on a RANS calculation. This can provide the mean flow, and also the decorrelation time through the typical time of turbulence dissipation. The tensor $\sigma$ can be expanded onto resolvent suboptimal modes, as in the present study.

## VI. CONCLUSION

In this paper we have proposed a stochastic modelling strategy of coherent structures in turbulent channel flows. By adding a stochastic non-linear forcing term, we obtain a refinement of the model proposed in Tissot et al. [41]. This forced model aims at improving consistency with the physical processes involved in these flows while maintaining the mathematical assumptions of the stochastic formulation. We used this model to explore the prediction abilities in the buffer and logarithmic layers in channels with friction Reynolds number equal to 180, 550 and 1000.

A central ingredient is the incorporation of a non-linear forcing representing coherent wave-wave interactions, which are essential in the self-sustaining processes of wall-bounded turbulence to generate streamwise vortices for large-scale structures in the logarithmic layer. In the model predictions, such forcing is central in the buffer layer, consistently with the Hamilton-Kim-Waleffe scenario [5]. Moreover we have shown that it is also central for large-scale structures in the logarithmic layer, which is in line with a large body of evidence in the literature supporting that self-sustaining processes are in action in the logarithmic layer. In addition to this, the model predicts a significant effect of incoherent turbulence on large log-layer structures, showing that it is crucial to model this effect in order to obtain good predictions. The model that we propose takes into account stochastic transport by this incoherent velocity field leading to an improvement compared to the state-of-the-art $\nu_t$-resolvent analysis. The present model requires a knowledge of the statistics of the turbulent field, and there are some possible open directions concerning its specification.

The modelling ingredients act on several physical mechanisms. First the stochastic diffusion induced by incoherent velocity field, unlike an eddy viscosity, has the shape of a full tensor with in particular $\langle u'v' \rangle$ off-diagonal components, which are defined consistently with the RMS profiles. A mean drift velocity takes into account the turbophoresis effect (effective transport from high to low turbulence regions), which is active in the buffer layer. A stochastic term representing the lift-up induced by the random incoherent velocity field is explicitly taken into account.

In addition, we have brought technical improvements by ensuring the decorrelation of the incoherent component with the solution through an iterative procedure, and we have proposed an efficient computation of FSLM by reformulating it as a singular value decomposition problem. With these effects taken into account in the model, we obtained an agreement between model predictions and turbulent fluctuations at various wall-normal positions.

In summary, the proposed model incorporates features of resolvent analysis, via the forcing resulting from non-linear wave-wave interactions, to the stochastic formalism introduced in our previous work [41], combining hence the benefits of the two approaches. The stochastic framework can be seen as a refined model of incoherent turbulence on large-scale structures, which not only includes the standard additional diffusion present in eddy-viscosity models [69, 70], but also involves all aspects of stochastic transport at the lengthscale of interest. The study shows that FSLM seems to carry advantageous features for reduced-order modelling of coherent structures in turbulent flows.
The authors acknowledge the support of the ERC EU project 856408 STUOD and the support of Rennes Métropole.

[1] S. J. Kline, W. C. Reynolds, F. A. Schraub, and P. W. Runstadler. The structure of turbulent boundary layers. *Journal of Fluid Mechanics*, 30(04):741-773, 1967.

[2] C. R. Smith and S. P. Metzler. The characteristics of low-speed streaks in the near-wall region of a turbulent boundary layer. *Journal of Fluid Mechanics*, 129:27-54, 1983. doi:10.1017/S0022112083000634.

[3] J. Jiménez and P. Moin. The minimal flow unit in near-wall turbulence. *Journal of Fluid Mechanics*, 225(213-240), 1991.

[4] R. L. Panton. Overview of the self-sustaining mechanisms of wall turbulence. *Progress in Aerospace Sciences*, 37(4): 341-383, 2001.

[5] J. M. Hamilton, J. Kim, and F. Waleffe. Regeneration mechanisms of near-wall turbulence structures. *Journal of Fluid Mechanics*, 287:317-348, 1995. doi:10.1017/S0022112095000978.

[6] T. Ellingsen and E. Palm. Stability of linear flow. *The Physics of Fluids*, 18(4):487-488, 1975. doi:10.1063/1.861156.

[7] L. Brandt. The lift-up effect: The linear mechanism behind transition and turbulence in shear flows. *European Journal of Mechanics-B/Fluids*, 47:80-96, 2014.

[8] O. Flores and J. Jiménez. Hierarchy of minimal flow units in the logarithmic layer. *Physics of Fluids*, 22(7):071704, 2010. doi:http://dx.doi.org/10.1063/1.3464157. URL http://scitation.aip.org/content/aip/journal/pof2/22/7/10.1063/1.3464157.

[9] A. J. Smits, B. J. McKeon, and I. Marusic. High-Reynolds number wall turbulence. *Annual Review of Fluid Mechanics*, 43(1):353-375, 2011. doi:10.1146/annurev-fluid-122109-160753.

[10] S. B. Pope. *Turbulent Flows*. Cambridge University Press, 2000. doi:10.1017/CBO9780511840531.

[11] Y. Hwang and C. Cossu. Self-sustained process at large scales in turbulent channel flow. *Physical Review Letters*, 105: 044505, Jul 2010. doi:10.1103/PhysRevLett.105.044505. URL https://link.aps.org/doi/10.1103/PhysRevLett.105.044505.

[12] C. Cossu and Y. Hwang. Self-sustaining processes at all scales in wall-bounded turbulent shear flows. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 375(2089):20160088, 2017.

[13] A. Lozano-Durán, H. J. Bae, and M. P. Encinar. Causality of energy-containing eddies in wall turbulence. *Journal of Fluid Mechanics*, 882:A2, 2020. doi:10.1017/jfm.2019.801.

[14] H. J. Bae, A. Lozano-Durán, and B. J. McKeon. Nonlinear mechanism of the self-sustaining process in the buffer and logarithmic layer of wall-bounded flows. *Journal of Fluid Mechanics*, 914:A3, 2021. doi:10.1017/jfm.2020.857.

[15] D. Barkley. Linear analysis of the cylinder wake mean flow. *Europhysics Letters (EPL)*, 75(5):750-756, sep 2006. doi: 10.1209/epl/i2006-10168-7.

[16] P. J. Schmid and D. S. Henningson. *Stability and transition in shear flows*, volume 142. Springer-Verlag, 2001.

[17] L. N. Trefethen and M. Embree. *Spectra and pseudospectra: the behavior of nonnormal matrices and operators*. Princeton University Press, 2005.

[18] M. R. Jovanović and B. Bamieh. Componentwise energy amplification in channel flows. *Journal of Fluid Mechanics*, 534: 145-183, 2005. doi:10.1017/S0022112005004295.

[19] B. J. McKeon and A. S. Sharma. A critical-layer framework for turbulent pipe flow. *Journal of Fluid Mechanics*, 658: 336-382, 2010.

[20] F. Gómez, A. S. Sharma, and H. M. Blackburn. Estimation of unsteady aerodynamic forces using pointwise velocity data. *Journal of Fluid Mechanics*, 804:R4, 2016. doi:10.1017/jfm.2016.546.

[21] S. Symon, D. Sipp, and B. J. McKeon. A tale of two airfoils: resolvent-based modelling of an oscillator versus an amplifier from an experimental mean. *Journal of Fluid Mechanics*, 881:51-83, 2019. doi:10.1017/jfm.2019.747.

[22] E. Martini, A. V. G. Cavalieri, P. Jordan, A. Towne, and L. Lesshafft. Resolvent-based optimal estimation of transitional and turbulent flows. *Journal of Fluid Mechanics*, 900:A2, 2020. doi:10.1017/jfm.2020.435.

[23] A. Towne, A. Lozano-Durán, and X. Yang. Resolvent-based estimation of space-time flow statistics. *Journal of Fluid Mechanics*, 883:A17, 2020. doi:10.1017/jfm.2019.854.

[24] F. R. Amaral, A. V. G. Cavalieri, E. Martini, P. Jordan, and A. Towne. Resolvent-based estimation of turbulent channel flow using wall measurements. *Journal of Fluid Mechanics*, 927:A17, 2021. doi:10.1017/jfm.2021.764.

[25] L. Franceschini, D. Sipp, and O. Marquet. Mean- and unsteady-flow reconstruction with one or two time-resolved measurements. *arXiv preprint arXiv:2102.03839*, 2021.

[26] C. Leclercq, F. Demourant, C. Poussot-Vassal, and D. Sipp. Linear iterative method for closed-loop control of quasiperiodic flows. *Journal of Fluid Mechanics*, 868:26-65, 2019. doi:10.1017/jfm.2019.112.

[27] L. I. Abreu, A. V. G. Cavalieri, P. Schlatter, R. Vinuesa, and D. S. Henningson. Resolvent modelling of near-wall coherent structures in turbulent channel flow. *International Journal of Heat and Fluid Flow*, 85:108662, 2020.

[28] S. Beneddine, D. Sipp, A. Arnault, J. Dandois, and L. Lesshafft. Conditions for validity of mean flow stability analysis. *Journal of Fluid Mechanics*, 798:485-504, 2016. doi:10.1017/jfm.2016.331. URL https://www.cambridge.org/core/article/conditions-for-validity-of-mean-flow-stability-analysis/2B42614C79F7EE9860AA56FB3DF1AB39.

[29] W. C. Reynolds and A. K. M. F. Hussain. The mechanics of an organized wave in turbulent shear flow. Part
3. Theoretical models and comparisons with experiments. *Journal of Fluid Mechanics*, 54(2):263-288, 1972. doi: 10.1017/S0022112072000679.

[30] R. D. Cess. A survey of the literature on heat transfer in turbulent tube flow. *Research Report No. 8-0529-R24*, 1958.

[31] Y. Hwang and C. Cossu. Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow. *Journal of Fluid Mechanics*, 664:51-73, 2010. doi:10.1017/S0022112010003629.

[32] P. Morra, O. Semeraro, D. S. Henningson, and C. Cossu. On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows. *Journal of Fluid Mechanics*, 867:969-984, 2019. doi:10.1017/jfm.2019.196.

[33] S. Symon, S. J. Illingworth, and I. Marusic. Large-scale structures predicted by linear models of wall-bounded turbulence. *Journal of Physics: Conference Series*, 1522(1):012006, 2020.

[34] P. Morra, P. A. S. Nogueira, A. V. G. Cavalieri, and D. S. Henningson. The colour of forcing statistics in resolvent analyses of turbulent channel flows. *Journal of Fluid Mechanics*, 907:A24, 2021. doi:10.1017/jfm.2020.802.

[35] S. Symon, S. J. Illingworth, and I. Marusic. Energy transfer in turbulent channel flows and implications for resolvent modelling. *Journal of Fluid Mechanics*, 911, 2021.

[36] S. Symon, A. Madhusudanan, S. J. Illingworth, and I. Marusic. On the use of eddy viscosity in resolvent analysis of turbulent channel flow. *arXiv*, 2022. doi:10.48550/ARXIV.2205.11216. URL https://arxiv.org/abs/2205.11216.

[37] P. A. S. Nogueira, P. Morra, E. Martini, A. V. G. Cavalieri, and D. S. Henningson. Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow. *Journal of Fluid Mechanics*, 908:A32, 2021. doi:10.1017/jfm.2020.918.

[38] V. Gupta, A. Madhusudanan, M. Wan, S. J. Illingworth, and M. P. Juniper. Linear-model-based estimation in wall turbulence: improved stochastic forcing and eddy viscosity terms. *Journal of Fluid Mechanics*, 925:A18, 2021. doi: 10.1017/jfm.2021.671.

[39] A. Zare, M. R. Jovanović, and T. T. Georgiou. Colour of turbulence. *Journal of Fluid Mechanics*, 812:636-680, 2017. doi:10.1017/jfm.2016.682.

[40] A. Zare, T. T. Georgiou, and M. R. Jovanović. Stochastic dynamical modeling of turbulent flows. *Annual Review of Control, Robotics, and Autonomous Systems*, 3, 2019.

[41] G. Tissot, A. V. G. Cavalieri, and E. Mémin. Stochastic linear modes in a turbulent channel flow. *Journal of Fluid Mechanics*, 912:A51, 2021. doi:10.1017/jfm.2020.1168.

[42] E. Mémin. Fluid flow dynamics under location uncertainty. *Geophysical & Astrophysical Fluid Dynamics*, 108(2):119-146, 2014.

[43] P. Chandramouli, D. Heitz, S. Laizet, and E. Mémin. Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty. *Computers & Fluids*, 168:170-189, 2018. ISSN 0045-7930. doi:https://doi.org/10.1016/j.compfluid.2018.04.001. URL http://www.sciencedirect.com/science/article/pii/S0045793018301890.

[44] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part I Random transport and general models. *Geophysical and Astrophysical Fluid Dynamics*, 111(3):149-176, April 2017. doi: 10.1080/03091929.2017.1310210. URL https://hal.inria.fr/hal-01391420.

[45] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading. *Geophysical and Astrophysical Fluid Dynamics*, 111(3):177-208, April 2017. doi: 10.1080/03091929.2017.1312101. URL https://hal.inria.fr/hal-01391476.

[46] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions. *Geophysical and Astrophysical Fluid Dynamics*, 111(3):209-227, April 2017. doi:10.1080/03091929.2017.1312102. URL https://hal.inria.fr/hal-01391484.

[47] B. Chapron, P. Dérian, E. Mémin, and V. Resseguier. Large scale flows under location uncertainty: a consistent stochastic framework. *Quarterly Journal of the Royal Meteorological Society*, 144(710):251-260, 2018. doi:10.1002/qj.3198. URL https://hal.inria.fr/hal-01629898.

[48] W. Bauer, P. Chandramouli, B. Chapron, L. Li, and E. Mémin. Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach. *Journal of Physical Oceanography*, February 2020. doi:10.1175/JPO-D-19-0164.1. URL https://hal.inria.fr/hal-02398521.

[49] W. Bauer, P. Chandramouli, L. Li, and E. Mémin. Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models. *Ocean Modelling*, 151:1-50, 2020. doi:10.1016/j.ocemod.2020.101646. URL https://hal.inria.fr/hal-02666147.

[50] B. Pinier, E. Mémin, S. Laizet, and R. Lewandowski. A stochastic flow approach to model the mean velocity profile of wall-bounded flows. *Physical Review E*, 2019. URL https://hal.inria.fr/hal-01947662.

[51] Y. Yang and E. Mémin. High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar. *Tellus A*, page 19, April 2017. URL https://hal.inria.fr/hal-01500140.

[52] Y. Yang and E. Mémin. Estimation of physical parameters under location uncertainty using an Ensemble²-Expectation-Maximization algorithm. *Quarterly journal of the royal meteorological society*, pages 1-47, November 2018. doi: 10.1002/qj.3438. URL https://hal.inria.fr/hal-01944730.

[53] P. Chandramouli, E. Mémin, and D. Heitz. 4D large scale variational data assimilation of a turbulent flow with a dynamics error model. *Journal of Computational Physics*, 412, July 2020. doi:10.1016/j.jcp.2020.109446. URL https://hal.inria.fr/hal-02547763.

[54] V. Resseguier, E. Mémin, D. Heitz, and B. Chapron. Stochastic modelling and diffusion modes for proper orthogonal decomposition models and small-scale flow analysis. *Journal of Fluid Mechanics*, 826:888-917, 2017.

[55] V. Resseguier, A. M. Picard, E. Mémin, and B. Chapron. Quantifying truncation-related uncertainties in unsteady fluid dynamics reduced order models. *SIAM/ASA Journal on Uncertainty Quantification*, 9(3):1152-1183, 2021.
596

[56] M. Reeks. The transport of discrete particles in inhomogeneous turbulence. *J. Aerosol Sci.*, 14(6):729-739, 1983.

597

[57] A. Towne, O. T Schmidt, and T. Colonius. Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. *Journal of Fluid Mechanics*, 847:821-867, 2018.

598

599

[58] E. Yim, P. Meliga, and F. Gallaire. Self-consistent triple decomposition of the turbulent flow over a backward-facing step under finite amplitude harmonic forcing. *Proceedings of the Royal Society A*, 475(2225):20190018, 2019.

600

601

[59] P. E. Plotter. Stochastic integration and differential equation. *Stochastic Modeling and Applied Probability*, 21, 2005.

602

[60] R. Moarref, A. S. Sharma, J. A. Tropp, and B. J. McKeon. Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels. *Journal of Fluid Mechanics*, 734:275-316, 2013. doi:10.1017/jfm.2013.457.

603

604

[61] J. H. M. Ribeiro, C.-A. Yeh, and K. Taira. Randomized resolvent analysis. *Physical Review Fluids*, 5:033902, Mar 2020. doi:10.1103/PhysRevFluids.5.033902. URL https://link.aps.org/doi/10.1103/PhysRevFluids.5.033902.

605

606

[62] E. Martini, D. Rodríguez, A. Towne, and A.V.G. Cavalieri. Efficient computation of global resolvent modes. *Journal of Fluid Mechanics*, 919:A3, 2021. doi:10.1017/jfm.2021.364.

607

608

609

[63] J. Jiménez. Near-wall turbulence. *Physics of Fluids (1994-present)*, 25(10):101302, 2013. doi: http://dx.doi.org/10.1063/1.4824988.

610

611

[64] S. Kadri Harouna and E. Mémin. Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling. *Computers and Fluids*, 156:456-469, August 2017. doi:10.1016/j.compfluid.2017.08.017. URL https://hal.inria.fr/hal-01394780.

612

613

[65] J. F. Gibson, F. Reetz, S. Azimi, A. Ferraro, T. Kreilos, H. Schrobsdorff, M. Farano, A. F. Yesil, S. S. Schütz, M. Culpo, and T. M. Schneider. Channelflow 2.0. *Manuscript in preparation*, 2019.

614

615

[66] A. V. G. Cavalieri, P. Jordan, and L. Lesshafft. Wave-packet models for jet dynamics and sound radiation. *Applied Mechanics Reviews*, 71(2), 2019.

616

617

[67] A. V. G. Cavalieri, D. Rodriguez, P. Jordan, T. Colonius, and Y. Gervais. Wavepackets in the velocity field of turbulent jets. *Journal of Fluid Mechanics*, 730:559-592, 9 2013. ISSN 1469-7645. doi:10.1017/jfm.2013.346. URL http://journals.cambridge.org/article_S0022112013003467.

618

619

620

[68] J. C. Del Álamo and J. Jiménez. Linear energy amplification in turbulent channels. *Journal of Fluid Mechanics*, 559: 205-213, 2006.

621

622

623

[69] E. Pickering, G. Rigas, O. T. Schmidt, D. Sipp, and T. Colonius. Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets. *Journal of Fluid Mechanics*, 917:A29, 2021. doi:10.1017/jfm.2021.232.

624

625

[70] C. Cossu, G. Pujals, and S. Depardon. Optimal transient growth and very large-scale structures in turbulent boundary layers. *Journal of Fluid Mechanics*, 619:79-94, 2009. doi:10.1017/S0022112008004370.