Gilles Tissot

† André

V G Cavalieri

Étienne Mémin

Stochastic linear modes in a turbulent channel flow

This study is focused on the prediction of coherent structures, propagating within a turbulent channel flow. We propose a derivation of the linearised problem based on a stochastic formulation of the Navier-Stokes equations. It consists in considering the transport of quantities by a resolved velocity (i.e. solution of the model) perturbed by a Brownian motion which models the unresolved turbulent fluctuations over the timeaveraged field, here thought of as the underlying background turbulence. The associated linearised model, considering the mean velocity profile as given, predicts linear solutions evolving within a corrected mean velocity field and perturbed by modelled background turbulence. Two ways to define the statistics of the Brownian motion are proposed and compared: one based on full simulation data, and the second, data-free, based on preliminary predictions from resolvent analysis. The technique is applied on turbulent channel flows at friction Reynolds number Re τ = 180 and Re τ = 550, and predictions are compared to direct numerical simulation results. We show that the principal components of an ensemble of solutions of this stochastic linearised system are able to represent the leading proper orthogonal decomposition (SPOD) modes with a similar accuracy than optimal responses coming form resolvent analysis with eddy viscosity model at scales where strong production occurs. For the other scales, receiving energy by non-linear redistribution, the present strategy improves the prediction. Moreover, the second mode is systematically well predicted over all scales. This behaviour is understood by the ability of the stochastic modelling to model positive and negative inter-scale energy transfers through stochastic diffusion and random stochastic transport, while the eddy viscosity term in resolvent analysis is purely diffusive.

Introduction

Waves evolving within turbulent flows, reminiscent of modes built from linear stability theory, are of particular interest in fluid flow modelling. These structured patterns are related to large coherent structures of velocity and pressure fluctuations, which can induce drag, mixing or acoustic emission. In the case of laminar flows, eigensolutions of the problem linearised over the steady state is an appropriate manner of dealing with infinitesimal perturbations growing in time and/or in space [START_REF] Pier | On the frequency selection of finite-amplitude vortex shedding in the cylinder wake[END_REF]). In the case of turbulent flows, even if the underlying hypotheses of infinitesimal disturbances are broken, analysis of similar linearised problems has been shown to be predictive under specific circumstances: linearising over the mean flow instead of the steady state [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF], and agreement is often obtained when there is a preferred amplification mechanism that dominates the flow response [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. Parallel free shear flows such as jets [START_REF] Cavalieri | Wavepackets in the velocity field of turbulent jets[END_REF] and mixing layers are subject to these properties, at least in the upstream regions exhibiting a strong dominance of the Kelvin-Helmholtz mode.

In recent years, resolvent analysis [START_REF] Schmid | Stability and transition in shear flows[END_REF][START_REF] Trefethen | Spectra and pseudospectra: the behavior of nonnormal matrices and operators[END_REF] has become extremely popular due to its ability to model the response of the linearised system at a given frequency to an unknown forcing [START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF][START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]. When dealing with turbulent flows, such forcing represents the otherwise omitted non-linear term. This technique was successfully able to represent coherent structures present in regions where linear stability fails. We can cite as examples wall-bounded turbulence [START_REF] Luhar | On the structure and origin of pressure fluctuations in wall turbulence: predictions based on the resolvent analysis[END_REF][START_REF] Sharma | Low-dimensional representations of exact coherent states of the Navier-Stokes equations from the resolvent model of wall turbulence[END_REF] or downstream regions of jets [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF][START_REF] Lesshafft | Resolvent-based modeling of coherent wave packets in a turbulent jet[END_REF]. Beyond flow analysis, it has been used in the context of reduced-order modelling, data assimilation and control (Gómez et al. 2016a,b;[START_REF] Leclercq | Linear iterative method for closed-loop control of quasiperiodic flows[END_REF][START_REF] Symon | A tale of two airfoils: resolvent-based modelling of an oscillator versus an amplifier from an experimental mean[END_REF]Martini et al. 2020a;[START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF]. Resolvent analysis provides a basis for the forcing and a basis of associated responses sorted by response gain. Unfortunately, it does not give information on the non-linearities present in the flow. Resolvent analysis leads to an exact agreement with flow statistics if and only if the forcing is a Gaussian white noise [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF][START_REF] Cavalieri | Wave-packet models for jet dynamics and sound radiation[END_REF]. In this case, response modes match the eigenfunctions of the cross-spectral density (CSD) of the flow, referred to as spectral proper orthogonal decomposition (SPOD) modes. In the temporal domain, covariance matrices can be predicted from time-decorrelated stochastic forcing covariance matrices by solving a Lyapunov equation [START_REF] Farrell | Stochastic forcing of the linearized Navier-Stokes equations[END_REF]. Links between temporal and frequency domain formulations can be found in [START_REF] Farrell | Generalized stability theory. part i: Autonomous operators[END_REF] and in [START_REF] Dergham | Stochastic dynamics and model reduction of amplifier flows: the backward facing step flow[END_REF] for amplifier flows.

Obviously, in turbulent flows, especially in shear flows where strong inhomogeneities occur, non-linearities are not white noise. In recent studies, some attempts to introduce a coloured forcing are performed. Strategies to colour an additive noise of temporal models using techniques from fluid mechanics modelling and control theory are reviewed in [START_REF] Zare | Stochastic dynamical modeling of turbulent flows[END_REF]. In [START_REF] Zare | Colour of turbulence[END_REF], the covariance of the forcing in the temporal domain is found by an optimisation problem in order to match the empirical covariance matrix under a Lyapunov equation constraint and maximising an entropy criterion. In the frequency domain, [START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF] proposed to use the resolvent operator in order to link the CSD of the forcing with the CSD of the response and thus reconstruct flow statistics from partially observed measurements. Inverse problems can be defined in order to identify forcing projections in order to obtain the right flow statistics, as in [START_REF] Moarref | A low-order decomposition of turbulent channel flow via resolvent analysis and convex optimization[END_REF]. In a different approach, intensive data-processing enables to extract from direct numerical simulations the flow statistics of the non-linear term, considered as an external forcing of the linearised Navier-Stokes equations in order to predict flow statistics [START_REF] Nogueira | Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow[END_REF][START_REF] Morra | The colour of forcing statistics in resolvent analyses of turbulent channel flows[END_REF]. The latter approach leads to exact agreement with flow data, but due to the use of the full time-series forcing information, it cannot be considered as predictive techniques, but rather a data analysis tool. It would therefore be important to understand how to best model the unknown forcing in a general setting, allowing a refinement of the predictions of resolvent analysis with a better approximation of these underlying non-linearities in the flow.

Another improvement of standard resolvent analysis has been to introduce an eddy viscosity in the resolvent operator [START_REF] Kaiser | Prediction of the flow response of a turbulent flame to acoustic pertubations based on mean flow resolvent analysis[END_REF][START_REF] Morra | On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows[END_REF]. The choice of linearisation is arbitrary, and thus an additional term can be added to represent the turbulent diffusion. In this case the time-averaged flow becomes a fixed point of the Reynolds-averaged Navier-Stokes (RANS) equations, which, although not strictly required in the linearisation procedure, may be thought as an advantage of eddy-viscosity models. On the other hand, it should be mentioned that such eddy-viscosity models are derived in order to model the Reynolds stresses in the mean-flow equation, and, despite their good performance in some cases, there is no guarantee that eddy-viscosity models designed in the RANS setting are appropriate to model generalised Reynolds stresses at non-zero frequencies and wavenumbers, i.e. considering Fourier transforms instead of time averaging. This caveat is supported by the analysis of [START_REF] Symon | Energy transfer in turbulent channel flows and implications for resolvent modelling[END_REF], showing that the eddy viscosity models well the non-linear transfers of exceeding energy at scales where strong production occurs, but is unable to take into account these transfers for scales receiving energy from the others since eddy viscosity is purely diffusive. This is consistent with the framework for which it has been defined in [START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments[END_REF]; eddy viscosity has been introduced in the context of a triple decomposition, with the warning that the representation by an eddy viscosity of the oscillation of the background Reynolds stress due to the passage of organised disturbances is a priori valid only for low frequencies and weak oscillations of large scales.

In the present study, we propose to adopt a different strategy. We start with a stochastic version of the Navier-Stokes equations under location uncertainty, as proposed by [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]. This model corresponds to the stochastic transport of mass and momentum by a resolved time-differentiable velocity perturbed by a Brownian motion representing an incoherent unresolved turbulent velocity field. Compared to the deterministic Navier-Stokes equations, it exhibits three additional features: i) a stochastic diffusion term coming from the effect of the stochastic transport, ii) the Brownian perturbation of the transport velocity and iii) a correction of the resolved transport velocity, called drift correction, caused by the inhomogeneity of variance of the unresolved velocities. This model has been successfully used for large eddy simulation (LES) and data assimilation (Resseguier et al. 2017a;[START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF][START_REF] Chandramouli | Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty[END_REF]. In this context the coloured Brownian motion represents the unresolved velocity field, the resolved velocity is the LES-filtered field and the stochastic diffusion can be interpreted as a turbulent eddy viscosity tensor. The drift velocity has been shown to be important in wall-bounded flows to capture profiles in the transitional buffer region where strong inhomogeneity occur [START_REF] Pinier | A stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF].

In our context, we consider a triple decomposition where the unresolved velocity is an incoherent contribution of the turbulent velocity fluctuation over the time-averaged field, and the resolved velocity is split into the time-averaged field † and a low-amplitude coherent structure, solution of a linearised model. In the model, we take into account the transport of the linear solution (a linear wave) by the turbulent fluctuation, and we assume that the linear wave does not affect the stochastic baseflow since the latter is modelled (and not computed). Hence, we do not consider transport of the turbulent fluctuation by the linear solution. By expressing this model in the Fourier domain, linearising it over the mean velocity, and considering normal modes, we provide a stochastic linear model. We are thus able to obtain an ensemble of solutions, whose eigenvalue decomposition leads to modes comparable to SPOD. By doing this, we do not explicitly take into account the effect of non-linearities. We instead consider a linear wave evolving within a turbulent stochastic field with known (or assumed) properties. The obtained linearised model is thus inhomogeneous, with coloured forcing provided by stochastic transport and the modelled Brownian motion. The direct application is to model coherent structures in turbulent flows by predicting SPOD modes of the direct numerical simulation (DNS), as will be pursued here; however, the present model may also be applied to study the evolution of small disturbances introduced to turbulent flows, as in [START_REF] Iyer | Identifying eigenmodes of averaged small-amplitude perturbations to turbulent channel flow[END_REF]. A clear advantage of the method is that the user-specified stochastic field covariances are expressed in velocity and not in forcing as it is the case in standard additive-noise techniques. Thus, this facilitate its definition and eventual intensive dataprocessing [START_REF] Nogueira | Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow[END_REF][START_REF] Morra | The colour of forcing statistics in resolvent analyses of turbulent channel flows[END_REF] or related inverse problems [START_REF] Zare | Colour of turbulence[END_REF]. The proposed formalism ensures consistency between energy brought by the stochastic noise and the induced stochastic diffusion in a unified framework.

The remainder of the manuscript is organised as follows. In section 2 we present turbulent channel flow configurations used to test the method. In section 3 we present spectral POD and resolvent analysis that are standard modelling procedures to respectively extract and predict waves evolving within turbulent flows. In section 4.1 we present the formulation of the stochastic modelling under location uncertainty and its application for incompressible flows. In section 4.2 a novel stochastic linear model under location uncertainty is proposed. The model involves statistics of the assumed Brownian motion, and in section 4.3 we discuss modelling strategies adopted here. Finally in section 5, we apply this new model to channel flow configurations. The manuscript is completed with conclusions in section 6.

Flow configuration

In this paper, we consider incompressible turbulent channel flow at Re τ = 180 and Re τ = 550, with Re τ the friction Reynolds number based on the friction velocity u τ and the channel half-height h. We have carried out DNS with periodic boundary conditions in streamwise and spanwise directions, with box dimensions of (L x = 4π, L y = 2, L z = 2π) for Re τ = 180 (following [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF]) and (L x = 2π, Ly = 2, L z = π) for Re τ = 550 (the minimal dimensions for accurate channel statistics, as identified by Lozano-Durán & Jiménez 2014). Simulations were carried out using Channelflow 2.0 (Gibson et al. 2019). The time steps of the database are respectively ∆t = 0.5 and ∆t = 0.1 in outer units (based on bulk velocity and channel half-height), with total simulation duration given respectively by T = 1000 and T = 300. In wall units, the time steps correspond respectively to ∆t + = 5 and ∆t + = 3. The same datasets were used by [START_REF] Morra | The colour of forcing statistics in resolvent analyses of turbulent channel flows[END_REF] for an analysis of non-linear terms and their role in resolvent analysis. Validation results are presented in the cited work, and will not be repeated here for conciseness. The present simulations are able to accurately reproduce mean flow and root-mean-square velocities of the simulations of Del [START_REF] Del Álamo | Spectra of the very large anisotropic scales in turbulent channels[END_REF].

In the following, all variables are expressed in wall units, based on friction velocity and kinematic viscosity. We consider Cartesian coordinates x = (x, y, z) denoting streamwise, wall-normal and spanwise directions, respectively. The flow variables, i.e. the velocity components (u, v, w) and the pressure p, are split into a time average and fluctuation q = q + q , with q = (u, v, w, p) T , q = (U (y), 0, 0, p) T and q = (u , v , w , p ) T . We define as well the velocity vector u = ū + u = (u, v, w) T . Due to the invariance in the streamwise (x) and spanwise (z) directions by translational invariance of the mean flow, and also due to the periodic boundary conditions used in x and z directions, we define the Fourier decomposition with the following convention

q (x, y, z, t) = 1 2π α,β ∞ -∞
qα,β,ω (y)e i(αx+βz-ωt) dω.

(2.1) α and β are respectively the streamwise and spanwise wavenumbers and ω the angular frequency. qα,β,ω (y) is the space-time Fourier coefficient.

Standard approaches for data analysis and linearised modelling

3.1. Spectral POD

The frequency domain form of the proper orthogonal decomposition (POD), referred to as SPOD, is an orthonormal basis of modes that i) oscillate at a given frequency, ii) are perfectly coherent, and iii) are decorrelated from each other. It constitutes a procedure to extract coherent structures from data evolving at a given frequency. A complete description is performed in [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. SPOD modes are eigenfunctions of the cross spectral density (CSD) tensor S, obtained as

S(α, y, y , β, ω) = 1 L x L z ∞ -∞ Lx 0 Lz 0 C(x -x , y, y , z -z , τ )e i(ωτ -α(x-x )-β(z-z )) dτ d(x -x )d(z -z ), (3.1)
where C(xx , y, y , zz , tt ) is the two-point space-time correlation tensor. Homogeneity in x, z and t allows us to express the dependence in each direction using Fourier modes, as shown in equation (3.1). Thus, for each (α, β, ω) combination, the eigendecomposition

U (n) = u (n) 1 , . . . , u (n) N , (3.3)
where n is the index related to the data block. A discrete Fourier transform (DFT) is performed on each block

Û(n) = û(n) ω1 , . . . , û(n) ω N , (3.4) with û(n) ω k = N j=1 w j u (n) j e i2π(k-1)(j-1)/N ∆t, (3.5)
and w j are weights of a window function designed to reduce effects of non-periodicity of the signal in the data block. The window is normalised such that N i=1 w 2 i = 1 to be consistent with the energy contained in the signal. Then, a POD is performed onto the ensemble of Fourier components at a given frequency such that

S ω k W Φ SPOD j,ω k = λ j,ω k Φ SPOD j,ω k , (3.6)
with the estimate of the cross-spectral density matrix

S ω k = 1 ∆tM M n=1 û(n) ω k û(n) ω k * , (3.7)
and W a matrix of weights defining the spatial inner-product (a, b) W = a * W b, numerical approximation of (a(y), b(y)) y , with • * denoting the transpose-conjugate operation. Functions Φ SPOD j,ω k are the SPOD modes and λ j,ω k are the associated eigenvalues. In our numerical tests, we use N = 256 (resp. N = 512) at Re τ = 180 (resp. Re τ = 550) with an overlap of 3 4 N and w j = 2 sin 2 (π(j -1)/N )). We can see that SPOD modes are principal components of an ensemble of solutions defined in the Fourier domain. This random variability between ensemble members û(n) ω of the same statistically stationary flow reflects the effect of turbulence on waves, which cannot be properly represented by a single deterministic Fourier mode. Indeed, a deterministic purely coherent wave not affected by the turbulence would result in identical values of û(n) ω and then in a single SPOD mode equal to the Fourier coefficient. We will see that the stochastic linear problem presented in section 4.2 shares common structures with the SPOD that makes them consistent and comparable.

Resolvent analysis

Resolvent analysis aims at exploring the response of a linearised system to external forcing, interpreted as unknown non-linearities when analysing turbulent flows [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF][START_REF] Hwang | Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow[END_REF]. We start with the incompressible Navier-Stokes equations

   ∂u ∂t + (u • ∇) u = -∇p + 1 Re ∆u ∇ • u = 0.
(3.8)

Expressing equation (3.8) in the Fourier domain and splitting the solution into the parallel mean flow U (y) and fluctuations u leads to the system

       A(•) + D(•) ∂U ∂y 0 iα 0 A(•) + D(•) 0 ∂• ∂y 0 0 A(•) + D(•) iβ iα ∂• ∂y iβ 0            ûα,β,ω vα,β,ω ŵα,β,ω pα,β,ω     =     fx fy fz 0     ,
(3.9)

with the advection operator A(•) = -iω + iαU , the diffusion operator

D(•) = -1 Re -α 2 + ∂ 2 • ∂y 2 -β 2 , Re the Reynolds number and fx = F α,β,ω u • ∇u -u • ∇u fy = F α,β,ω u • ∇v -u • ∇v fz = F α,β,ω u • ∇w -u • ∇w , (3.10)
where F α,β,ω (•) denotes the space-time Fourier transform. Equation (3.9) can be written more compactly as (3.11) where A α,β,ω is the linearised Navier-Stokes operator under the considered assumptions, and E is a diagonal matrix with 1 and 0 values for velocity and pressure components respectively. L -1 α,β,ω = (A α,β,ω -iωE) -1 is the resolvent of the operator A α,β,ω . Formally, the non-linear term N α,β,ω (u ) can be considered as a right-hand side of equations (3.11). However, the Fourier transform of a product is a convolution, which results in an integral over all frequency-wavenumbers [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF][START_REF] Cavalieri | Wave-packet models for jet dynamics and sound radiation[END_REF], which is incompatible with a prediction for an "isolated" frequency-wavenumber combination; the non-linear term couples the various frequencies and wavenumbers in the Navier-Stokes system in triadic interactions. The goal of resolvent analysis is to explore the response of the linear operator L α,β,ω to non-linearities treated as an external forcing. In what follows we will consider the discretised form of equation (3.11), which, with the addition of an external forcing, becomes L α,β,ω qα,β,ω = B fα,β,ω .

(A α,β,ω -iωE) qα,β,ω = L α,β,ω qα,β,ω = N α,β,ω (u ),
(3.12)

Since non-linear terms act only on the momentum equations in (3.9), and since SPOD modes can be calculated separately once the velocity is known, we define input and output matrices respectively as

B =     I 0 0 0 I 0 0 0 I 0 0 0     and H =   I 0 0 0 0 I 0 0 0 0 I 0   .
(3.13)

The output matrix H thus extracts only the velocity components, which may be combined in an inner product whose norm is the kinetic energy. The response to a harmonic forcing fα,β,ω is H qα,β,ω = HL -1 α,β,ω B fα,β,ω . The frequency response operator is thus defined from the resolvent operator, L -1 α,β,ω , as R α,β,ω = HL -1 α,β,ω B. We then search for the forcing that maximises max fα,β,ω

H qα,β,ω 2 W fα,β,ω 2 W = R α,β,ω fα,β,ω 2 W fα,β,ω 2 W . (3.14)
The maximisation of the Rayleigh quotient (3.14) can be achieved by means of the following singular value decomposition,

W 1/2 R α,β,ω W -1/2 = U r Σ r V * r , (3.15)
where W 1/2 is defined by the Cholesky decomposition W = W 1/2 (W 1/2 ) * with W defined in section 3.1, U r = (U r,1 , . . . , U r,N ) and V r = (V r,1 , . . . , V r,N ) are orthonormal matrices and Σ r = diag(s 1 , . . . , s N ) is a diagonal matrix. We define the optimal forcing modes as,

Ψ resolvent i = W -1/2 V r,i
, and the associated optimal response modes (also referred to as resolvent modes in the following) as, Φ resolvent i = W -1/2 U r,i . The singular values, s i , diagonal elements of Σ r sorted in descending order, indicate the gains associated with forcing-response mode pairs.

As shown in [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF], resolvent modes and SPOD correspond if and only if the forcing in the data correspond to a Gaussian white noise. However, if there is a separation in gain between the optimal forcing-response pair and suboptimal ones, even coloured forcing leads to a dominance of the optimal response in the flow CSD, unless the forcing is strongly biased towards suboptimals [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF][START_REF] Cavalieri | Wave-packet models for jet dynamics and sound radiation[END_REF]. A way to account for a part of this forcing is to consider a triple decomposition of the velocity field in the line of the work of [START_REF] Reynolds | Stability of turbulent channel flow, with application to Malkus's theory[END_REF]; [START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments[END_REF] where the Reynolds stresses induced by the incoherent part of the turbulent velocity through ensemble averaging is modelled by an eddy-viscosity in the linearised operator, noted L νt α,β,ω . It corresponds to a modification of the system (3.8), where the diffusion operator 1

Re ∆u is replaced by ∇

•(( 1 Re +ν t (y)) ∇u + ∇u T ), leading to the system         A(•) + D νt (•) ∂U ∂y -iα ∂ν ∂y 0 iα 0 A(•) + D νt (•) -∂ν ∂y ∂• ∂y 0 ∂• ∂y 0 -iβ ∂ν ∂y A(•) + D νt (•) iβ iα ∂• ∂y iβ 0             ûα,β,ω vα,β,ω ŵα,β,ω pα,β,ω     =     fx fy fz 0     , (3.16) with D νt (•) = -1 Re + ν t -α 2 + ∂ 2 • ∂y 2 -β 2 -∂ν t ∂y ∂• ∂y .
In this work, we choose to fix the turbulent viscosity parameter, ν t , with the Cess's model [START_REF] Cess | A survey of the literature on heat transfer in turbulent tubeflow[END_REF] as in Del Álamo & [START_REF] Del Álamo | Linear energy amplification in turbulent channels[END_REF] and [START_REF] Pinier | A stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF]:

ν t = 1 Re 1 2 1 + κ 2 Re 2 τ 9 (1 -y 2 ) 2 (1 + 2y 2 ) 2 (1 -e -y + A ) 2 1 2 - 1 2 , ( 3.17)
where y + = Re τ (1 -|y|), κ = 0.426 is the von Kármán constant and A = 25.5 is a constant chosen consistently following [START_REF] Pinier | A stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF]. The corresponding resolvent analysis will be referred to as eddy-viscosity resolvent analysis. For compactness and to avoid confusions, standard resolvent analysis will be noted ν-resolvent and eddy-viscosity resolvent will be noted ν t -resolvent. Recent results have shown that the inclusion of an eddy viscosity in the linearised operator improves the agreement between resolvent modes and simulation data [START_REF] Illingworth | Estimating large-scale structures in wall turbulence using linear models[END_REF][START_REF] Morra | On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows[END_REF], as the eddy viscosity models part of the forcing statistics [START_REF] Morra | The colour of forcing statistics in resolvent analyses of turbulent channel flows[END_REF].

In any case, the resolvent modes constitute a basis where any forcing can be projected onto the forcing modes and thus we can deduce the response to this forcing. Extracting statistics of the non-linear term is a hard task in itself, and examples are shown by [START_REF] Zare | Colour of turbulence[END_REF], [START_REF] Nogueira | Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow[END_REF], [START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF] and [START_REF] Morra | The colour of forcing statistics in resolvent analyses of turbulent channel flows[END_REF], where DNS data was used for that matter. To avoid the need of detailed data for predictions, some modelling of the unknown forcing is necessary. We propose in this paper to extend such an analysis to a framework introducing a statistical modelling of the fluctuating velocities and a stochastic transport equation, with modelling of unresolved velocity as a Brownian motion. This is described next.

A stochastic linearised model for coherent structures in turbulent flows

Stochastic Navier-Stokes equations under location uncertainty

The starting point of our modelling approach is the stochastic model proposed by [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]. It relies on the decomposition of the Lagrangian transport of a particle by a resolved velocity u perturbed by a function of a Brownian motion B t , written as a stochastic Itô process. The displacement X(x, t) of a particle from time t 0 to t is thus written as

X(x, t) = X(x, t 0 ) + t t0 u(x, t)dt + t t0 σdB t , (4.1)
or in a more usual and compact differential form as

dX(x, t) = u(x, t)dt + σdB t . (4.2)
The term σdB t hides a spatial convolution in the domain Ω between the modelled bounded deterministic correlation kernel σ(x, x , t) and dB t the increment of the Brownian motion B t (i.e. dB t = B t+dt -B t ) such that its i th component at the spatial location x is given by the following kernel integral expression

(σdB t ) i x = Ω σ ij (x, x , t)dB j t (x ) dx . (4.3)
This formulation models the transport of physical quantities by a smooth resolved velocity (i.e. differentiable in time), perturbed by a noise that is assumed to be correlated and smooth in space (through the correlation kernel σ) but decorrelated and nondifferentiable in time. This ground assumption in Itô calculus is as well classically performed in standard stochastic forcing models, and is relaxed in [START_REF] Zare | Colour of turbulence[END_REF] for a time-domain formulation. The strong hypothesis of time-decorrelated multiplicative noise is clearly a limitation of the proposed formalism. However, introducing a time memory would require considering high-order terms in the stochastic Taylor expansions in Itô calculus [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] and is beyond the scope of the present paper.

The statistical behaviour of the noise is an information that has to be specified. The operator σ is built in order to respect the given covariance tensor of the unresolved part

E (σdB t ) i x (σdB t ) j x = Ω Ω σ ik (x, y, t)E dB k t (y)dB l t (y ) σ lj (x , y , t ) dy dy = Ω Ω σ ik (x, y, t) (δ(y -y )δ(t -t )δ kl ) dt σ lj (x , y , t ) dy dy = Ω σ ik (x, y, t)σ kj (x , y, t ) dy δ(t -t )dt ! = Q(x, x , t) δ(t -t )dt. (4.4)
with dB i t the i th component of dB t and E the expectation operator. We define as well the variance tensor a(x, t) = Q(x, x, t). (4.5)

We can remark that the tensor a(x, t) has the dimension of a diffusion. It corresponds to the so-called quadratic variation of the stochastic process.

The stochastic transport equation of a conserved quantity θ can be obtained through a stochastic version of the Reynolds transport theorem and the Itô-Wentzell formula [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]Resseguier et al. 2017a):

d t θ + ∇ • (θu d )dt + ∇ • (θσdB t ) = ∇ • 1 2 a∇θ dt, (4.6)
with d t θ the temporal increment of θ and the drift velocity

u d = u - 1 2 ∇ • a + σ(∇ • σ). (4.7)
The stochastic transport displacement described by equation

(4.6) is u dt + (-1 2 ∇ • a + σ(∇ • σ)) dt + σdB t .
It contains a non-differentiable part σdB t which is directly linked to the hypothesis (4.2). Moreover, there is an additional time-differentiable term -1 2 ∇ • a dt, called corrective drift velocity, which is active when the random fluctuations are inhomogeneous in space. This term corresponds to a turbophoresis correction introduced in the literature to represent laden particles drift due to turbulence inhomogeneity -see Resseguier et al. (2017d ) and references therein. The second corrective drift velocity term σ(∇ • σ) takes into account divergence of the noise, and vanishes in the incompressible case since by mass conservation we obtain ∇ • σ = 0 as a constraint on the Brownian terms. Additionally, the stochastic diffusion term ∇ • 1 2 a∇u , which is akin to a subgrid diffusion tensor associated to a generalized (matrix form) Boussinesq turbulent viscosity assumption, ensures energy conservation of the stochastic transport operator (Resseguier et al. 2017a). It represents the spreading of the expectation of the solution due to the randomness. The fundamental difference between standard eddy viscosity models and the stochastic diffusion is that eddy viscosity is designed in order to reproduce the mean flow in RANS equations based on an analogy with the molecular diffusion (so-called the Boussinesq assumption), while the stochastic diffusion is defined here consistently with the decomposition operated with no supplementary assumption.

Considering density and momentum for θ, with a constant density, reduces to a stochastic version of the Navier-Stokes equations for the resolved velocity u,

                 d t u + (u d • ∇) u dt + (σdB t • ∇) u = -∇ (p t dt + dp t ) + 1 Re ∇ • (∇u) dt + ∇ • 1 2 a∇u dt + 1 Re ∇ • (∇σdB t ) ∇ • u d = 0; ∇ • σ = 0, u d = u - 1 2 ∇.a.
(4.8) This stochastic system has been successfully used in several different applicative contexts such as reduced-order models (Resseguier et al. 2017d ), geophysical flow modelling and analysis (Resseguier et al. 2017a,b,c;[START_REF] Chapron | Large scale flows under location uncertainty: a consistent stochastic framework[END_REF]Bauer et al. 2020a,b) or for data assimilation techniques of coarse scale flow dynamics [START_REF] Chandramouli | Fast 3D flow reconstructions from 2D cross-plane observations[END_REF][START_REF] Chandramouli | 4D large scale variational data assimilation of a turbulent flow with a dynamics error model[END_REF] or large eddy simulation [START_REF] Kadri Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF]. Details to derive this system are reported in appendix A. Besides the fact that they are written in a differential form, equations (4.8) have noticeable differences compared to the deterministic Navier-Stokes equations. First, the transport displacement is formulated through the stochastic transport described previously, which incorporates a modified advection term driven by inhomogeneity of the random fluctuations, a multiplicative random term advecting the resolved component by the random fluctuations and a subgrid stochastic diffusion related to the fluctuation variance tensor. Secondly, there is an additional pressure term dp t that is the random part of the pressure term. This term is necessary in order to balance the martingale part of the equations, and may be understood as a random normal stress induced by the divergence-free condition. There is as well a contribution of the viscous forces induced by the unresolved velocity. Finally, the continuity condition is modified because applying mass conservation shows that the drift velocity u d is solenoidal, instead of u. We can note that ∇ • u = 1 2 ∇ • ∇ • a and vanishes only for some specific definitions of the noise (such as spatially homogeneous noise) and not in the general case. The corrective drift term is caused by turbulence inhomogeneity and has no reason to be solenoidal. For a zero noise, the deterministic solenoidal condition is recovered. The corrective drift term has been shown to play a fundamental role in the buffer region of turbulent boundary flows [START_REF] Pinier | A stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF], and it has been also demonstrated to allow the triggering of secondary circulations in geophysical flows (Bauer et al. 2020a).

It is important to remark that even if σdB t has a Gaussian probability distribution, it is involved in a multiplicative relation with the gradient of the solution, and thus this Gaussian property is not inherited by u. This is a fundamental difference with a strategy based on the coloration of an additive Gaussian noise.

Stochastic linearised model

The strategy proposed in this paper to model coherent structures propagating within a turbulent flow is to express the stochastic model (4.8) in the Fourier domain with a normal-mode Ansatz, to linearise it and to solve the obtained linear model subject to an ensemble of realisations of the Brownian motion. Such coherent structures are naturally modelled as waves by appealing to a Fourier decomposition of the field.

The normal-mode Ansatz assumes homogeneity in time (t) and in the streamwise (x) and spanwise (z) directions, i.e. all statistical quantities are invariant in these directions. Thus, we can write the mean flow as ū = (U (y) 0 0) T , the variance tensor as a(y) and the noise diffusion operator as σ(xx , y, y , zz , tt ).

In order to express the stochastic equations in the Fourier domain, we introduce the power spectral density of the stochastic process σdB t , defined as the Fourier transform of the (two-point two-time) covariance tensor over a small increment of time

∆ t+∆ t E (S σ ω ) ds = ∞ -∞ E t+∆ t (σdB s ) i x t +∆ t (σdB s ) j x e iω(t-t ) d(t -t ) = ∞ -∞ t+∆ t Q(x, x , s)ds δ(t -t )e iω(t-t ) d(t -t ) = t+∆ t Q(x, x , s)ds, (4.9)
with Q(x, x , t) defined in (4.4). We can remark that E (S σ ω ) is independent of ω. In other words, it is white noise with a flat spectral content in frequency. We define then the random variables dξ ω representative of the time Fourier transform of σdB t , which shares the same power spectral density:

dξ ω = Ω σ(x, x , t)dη ω (x ) dx , (4.10)
where dη ω is the time Fourier component of dB t . Thus, dξ ω is homogeneous to the Fourier component of a displacement with a covariance consistent with σdB t since the following relation holds.

E dξ i ω (x)dξ j ω (x ) = E (σdB t ) i x (σdB t ) j x = Q ij (x, x )dt. (4.11)
Then, we define dξ α,β,ω , the space-time Fourier transform of σdB t as

dξ α,β,ω = 1 -1
σα,β (y, y ) dη α,β,ω (y )dy , (4.12)

with σα,β (y, y ) the space Fourier transform of σ(x-x , y, y , z -z , t-t ) and dη α,β,ω the space-time Fourier transform of dB t . The corresponding highly oscillating velocity will be written formally ξα,β,ω = The linearisation of the stochastic Navier-Stokes equations leads to

         d t u + ( ūd • ∇) u dt + (u • ∇) ū dt + (σdB t • ∇) ū + (σdB t • ∇) u = -∇ (p t dt + dp t ) + 1 Re ∇ • (∇u ) dt + ∇ • 1 2 a∇u dt ∇ • u = 0; ∇ • σ = 0.
(4.13) In the linearisation, non-linear terms in u are neglected assuming that the wave has a small amplitude, even if it is carried by a turbulent flow. As usually stated in linear stability analysis, the mean velocity is assumed to be solution of the steady stochastic equations. For the steady solution, by vanishing the time variation, the "dt" terms and "dB t " terms can be separated in two equations. The former is similar to the deterministic case, while the stochastic part states balance between advection by σdB t and the random pressure term ensuring stationarity of the solution. Finally, the random viscous term 1 Re ∇ • (∇σdB t ), represents the viscous stress induced by the Brownian motion. Here it is assumed to not act on the linear wave, but only on the turbulent field, for which an implicit evolution equation is replaced by the stochastic modelling. Thus, the random viscous term is neglected in equation (4.13).

For further developments, we consider a spectral representation of the noise and expand the differential of the Brownian motion onto a basis such that

σdB t = α,β,j Φ B,α,β,j (y)e i(αx+βz) dζ j,t (4.14)
where the profiles Φ B,α,β,j (y) are eigenfunctions of σα,β , and ζ j,t are scalar Brownian motions. This allows us to introduce a stochastic extension of the classical normal-mode Ansatz of the form:

u (x, y, z, t) = ûα,β,ω (y)e i(αx+βz-ωt-j ω j ζj,t) σdB t = dξ α,β,ω (y)e i(αx+βz-ωt-j ω j ζj,t) (4.15)

and we set θ = αx+βz -ωtj ω j ζ j,t . In this stochastic Ansatz the wave phase is now random. Upon applying this stochastic normal-mode Ansatz to (4.13) and considering the Brownian terms resulting from time differentiation leads to a relationship on the Brownian terms proportional to e iθ dζ j,t (simplified in the following relation) such as

j iω j ûα,β,ω dζ j,t = F α,β,ω (σdB t • ∇u ) + ∇q j dζ j,t , (4.16)
with ∇ = (iα ∂ • /∂y iβ) T and q j resulting from the splitting of the random pressure term dp t = j q j (y)e iθ dξ j,t + dp r,t . Equation (4.16) states that the transport of the wave by the random velocity induces phase randomness and random pressure perturbations. The term F α,β,ω (σdB t • ∇u ), Fourier transform of the transport of the wave by σdB t , is a convolution over all frequency-wavenumbers. It is a consequence of the bi-linear structure of the equations with respect to (u , σdB t ). In order to reconstruct the temporal behaviour of the wave, this term should be explicitly computed by considering a frequency-wavenumber discretisation; or modelled by a random variable. If we consider only predictions of ûα,β,ω and of cross-spectral density tensors, and if we assume a time and scale separation between Brownian motion and coherent wave ûα,β,ω , thus the terms involved in (4.16) become decoupled from all the others and it is not necessary to use equation (4.16). In the following of the paper, we adopt these assumptions and we obtain the system

                     -iωû α,β,ω + iαU d ûα,β,ω + vα,β,ω ∂U ∂y + iαp α,β,ω + D(û α,β,ω ) = -( ξα,β,ω ) y ∂U ∂y -iωv α,β,ω + iαU d vα,β,ω + ∂ pα,β,ω ∂y + D(v α,β,ω ) = 0 -iω ŵα,β,ω + iαU d ŵα,β,ω + iβ pα,β,ω + D( ŵα,β,ω ) = 0 iαû α,β,ω + ∂v α,
D(•) = - 1 Re -α 2 + ∂ 2 • ∂y 2 -β 2 - 1 2 -α 2 a xx + iαa xy ∂• ∂y -αβa xz +iα ∂a yx • ∂y + ∂ ∂y a yy ∂• ∂y + iβ ∂a yz • ∂y -αβa zx + iβa zy ∂• ∂y -β 2 a zz .
(4.18)

Rearranging in matrix form leads to

       A(•) + D(•) ∂U ∂y 0 iα 0 A(•) + D(•) 0 ∂• ∂y 0 0 A(•) + D(•) iβ iα ∂• ∂y iβ 0           ûα,β,ω vα,β,ω ŵα,β,ω pα,β,ω    =    fx 0 0 0    , (4.19)
with the modified advection operator A(•) = -iω + iαU d , and f x = -( ξα,β,ω ) y ∂U ∂y .

The variable pα,β,ω is the space-time Fourier component of both random pressure dp r,t and time-differentiable pressure p. There is no distinction between them in the Fourier space, and thus they may be combined into a single unknown function. In these equations, ξα,β,ω is a random variable, and as a consequence the solution is a random variable as well. Equation (4.17) is a linear system to solve. For any triplet (ω, α, β), a probability distribution of solutions can be found, which could be interpreted by the fact that due to randomness, there is a probability for a linear wave to be sustained. This trait is absent in a deterministic linear stability problem since a solution has to respect the dispersion relationship.

Compared to the system (3.9), the forcing acts only on the streamwise component due to the parallel mean flow (hypothesis explored in [START_REF] Nogueira | Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow[END_REF], a stochastic diffusion is added and transport velocities are modified by corrective drift and stochastic contribution. Formally, all these modifications could be interpreted after rearrangement as a stochastic right-hand side of system (3.9). However, this forcing provides intrinsically a non-Gaussian forcing with space-time structure (through the stochastic transport) that we believe to be physically more relevant than an additive ad hoc (often Gaussian) random forcing term. The sparse nature of the forcing is due to the fact that only linear terms are conserved †, thus limiting the possible interactions.

Considering an ensemble of N ens realisations of Gaussian white noise ηα,β,ω in equation (4.17) leads to an ensemble of solutions, whose statistics, such as cross-spectral density, may be computed. As a final step, SPOD is performed on such an ensemble of solutions ûα,β,ω , leading to an orthonormal basis of stochastic linearised modes (SLM), with corresponding eigenvalues indicating the variance of each mode in the statistics predicted by the model.

It is important to outline, that compared to standard linear stability or resolvent analysis, the present strategy does not require an eigenvalue decomposition or the inversion of the resolvent operator. Instead, an ensemble of N ens linear system solutions are required, which are individually computationally less expensive. The computational cost complexity is thus proportional to N ens . This can be viewed as a clear advantage in terms of CPU and memory costs for large scale sparse systems. Moreover, since the left-hand side of equation (4.19) is the same for all ensemble members, cost reduction can be obtained using for instance block Krylov methods. Similarly as in resolvent analysis where efficient numerical procedures are proposed [START_REF] Monokrousos | Global threedimensional optimal disturbances in the blasius boundary-layer flow using time-steppers[END_REF][START_REF] Moarref | Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels[END_REF][START_REF] Brynjell-Rahkola | Computing optimal forcing using laplace preconditioning[END_REF][START_REF] Ribeiro | Randomized resolvent analysis[END_REF]Martini et al. 2020b), it is believed that clever optimisations could be considered in the future for SLM.

The underlying modelling assumption here is that a single wave ûα,β,ω , with low amplitude compared to the overall background turbulence, is assumed to be transported by the mean flow U (y) (as in deterministic linearised models) but also by a background turbulence, modelled here as Brownian motion. The influence of turbulence on the coherent wave is manifest in the drift velocity, in the additional diffusion associated with the derivatives of the tensor a(y) (akin to an eddy viscosity and in equilibrium with the energy brought by the noise), and in the forcing term in the x-momentum equation, where background turbulence forces the coherent response.

Choice of σ

In the stochastic model, the statistics of the unresolved velocity have to be defined through σ. It is through this tensor that information of the turbulent velocity field affecting the wave is incorporated in the modelling. We propose in this section two different ways of choosing σ depending whether data are available or not. The first is a best-case situation where SPOD modes of the velocity field are available. The second consists in approximating them based on resolvent analysis in a model-based perspective.

Based on SPOD

The present method is a best-case scenario where SPOD are used to define σ. This is of course cyclic since the predicted wave will be compared to SPOD as well. The nonlinear problem is not solved by the stochastic linear problem, and therefore we test here if provision of the right (at least as right as possible) statistics of the noise to the linear stochastic system leads to a good representation of the waves associated to the original non-linear system.

SPOD is a data-driven method that extracts an orthonormal basis for the velocity field at a given frequency whose modes are optimal in terms of energy content. It is thus the optimal basis with respect to a kinetic energy criterion to represent the turbulent fluctuations at a given frequency. Within this framework, the variance tensor a is chosen to match the CSD of the turbulent velocity field. Thus, an expansion onto SPOD modes scaled by the SPOD eigenvalues is performed

a(y) = NSPOD j=1 Φ SPOD j,αβω (y)λ SPOD j,αβω Φ SPOD j,αβω (y) * τ (4.20)
and the spectral representation of the noise in the form of (4.14) is given as ξα,β,ω = NSPOD j=1 Φ SPOD j,αβω (y) λ SPOD j,αβω ηj , (4.21)

where ηj is a centered standard white noise variable †, while τ is a decorrelation time consistent with dt, the characteristic decorrelation time constant of fluctuation components. Indeed, a has the dimension of a diffusion ([L 2 /T ]) and a stochastic diffusion time scale has to be set. Since we are focusing on single frequency solutions, we can consider that below a given time, the faster fluctuations are fast enough to be considered as decorrelated. We here arbitrarily choose such time proportional to the half wave period,

τ = K π ω . (4.22)
The constant K = 100 allows us to obtain an order of magnitude of a consistent with Cess's eddy viscosity. This value, selected by inspection of figure 1(a), has been determined in terms of order of magnitude and not based on a fine iterative tuning. A very low order of magnitude for K leads to disappearance of the effects of stochastic diffusion and drift velocity, and a very large order of magnitude leads to unphysical results dominated by the noise. The fact that SPOD modes are purely coherent modes decorrelated from each other, justifies the expansion (4.21) by randomising the expansion coefficients.

It should be remarked that this definition of the operator σ requires a complete timeresolved dataset in order to be able to compute the SPOD, which is per se an intensive post-processing procedure.

Based on resolvent analysis

As a second approach, we use ν-resolvent analysis to build the background turbulence statistics σ at the considered frequency-wavenumber combination. This amounts to considering an imperfect prediction based on the resolvent, with the simplifying assumption of non-linear terms as spatial white-noise, to build a first estimate of the turbulence statistics, which are subsequently taken to define the coloured Brownian motion statistics σ in our model. In that perspective, a distinction is made between the linear response to generalised Reynolds stresses modelling the background turbulence σdB t and the linear wave sustained by this background turbulence qα,β,ω .

With this underlying idea, we define, similarly as for the SPOD noise

a resolvent (y) =   A r N resolvent j=1 Φ resolvent j,αβω (y) s resolvent j,αβω 2 Φ resolvent j,αβω (y) * + A b I   τ (4.23) and ( ξα,β,ω ) resolvent = A 1 2 r N resolvent j=1 Φ resolvent j,αβω (y)s resolvent j,αβω ηj + A 1 2 b ηbackground . (4.24)
The constant A r is a global amplitude that is a free parameter. For reasons of consistency with the data-driven noise defined from SPOD, we choose A 1 2 r = λ SPOD 1,αβω /s resolvent 1,αβω . This sets the noise amplitude from data, but it could be considered as well as a single free parameter in a fully model-based perspective. The resolvent predictions were obtained without inclusion of an eddy viscosity to avoid a circular nature of the model; the resolvent † Numerically generated by the Matlab/Octave function randn(). only provides an estimate of the stochastic field, which will subsequently be used to force the waves and modify the linearised operator via the drift velocity and the stochastic diffusion. Moreover, it has been observed that an additional constant background noise is necessary to represent fluctuations in the middle of the channel (structures exemplified by wave 2 as defined in section 5), present with SPOD noise, but missing with resolvent noise. The constant A b = 1/τ is chosen in order to obtain comparable levels than Cess's model in the middle of the channel. The noise term ηbackground denotes here a centered standard white noise field. The modelling is designed to mimic the behaviour of SPOD noise, and refinement of the definition of σ is an open research direction.

Application to turbulent channel flows

Effect of the variance tensor

As a first analysis, we show here the effect of the tensor a, through the stochastic diffusion and the drift velocity. In figure 1(a), we compare 1 2 trace(a) 3

of wave 1 at Re τ = 550 with Cess' eddy viscosity ν t . Compared to Cess' model that leads to an eddy viscosity that increases monotonically with wall distance, so as to lead to the mean velocity profile of the channel, SPOD and resolvent noise lead to large diffusion in the buffer layer around y + ≈ 10, where the SPOD modes peak. It was observed that the eddy viscosity/stochastic diffusion in the channel center plays a significant role in the SLM mode's shape. Indeed, SPOD noise leads to good results, but resolvent noise without adding the background homogeneous noise showed a faster decay in the channel center, which we believe to be the reason of failure. This is the reason why we have added an additional homogeneous noise in the resolvent noise equations (4.23) and (4.24), such that the level of induced stochastic diffusion is comparable to ν t . The effect of this noise can be shown for y + > 100 in the resolvent noise, and that constitutes the main difference between both noises. Let us note however that the stochastic model should not be interpreted as more diffusive than Cess' model. As a matter of fact, the noise which acts globally as a backscatering term is counterbalanced by the (negative) stochastic diffusion energy.

In figure 1(b), we compare the mean flow profile U (y) with the mean drift velocity profile U d (y) at Re τ = 550 with the SPOD noise of wave 1. The main effect of the corrective drift is around y + = 10, which is consistent with the study of [START_REF] Pinier | A stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF] showing that the stochastic transport formalism has a large influence in the buffer layer.

Comparison with SPOD modes from the DNS databases

The linear stochastic modes are compared with results of resolvent analysis, described in section 3.2. These models are discretised with 128 (resp. 256) Chebyshev polynomials at Re τ = 180 (resp. Re τ = 550), the same number of polynomials used in the DNS. The accuracy of SLM in modelling coherent structures in turbulent channel flows is evaluated by comparing results to SPOD modes from the two channels. The two choices of noise presented in section 4.3 are tested. For all the results, the two leading modes from SPOD, resolvent analysis and the linearised stochastic model are considered. They are referred to respectively as mode 1 and mode 2. Higher order modes are not shown in order to avoid any misinterpretation due to numerical convergence of the SPOD, which requires very long time series for higher order modes [START_REF] Lesshafft | Resolvent-based modeling of coherent wave packets in a turbulent jet[END_REF]. Stochastic linearised modes are obtained using an ensemble of N ens = 10000 members. In practice, 1000 members are enough to obtain robust results, and we have increased this value to eliminate any doubt of convergence.

Figures 2 and3 show mode amplitudes for Re τ = 180 and Re τ = 550 respectively. Recall that the wall-normal coordinate y is expressed in inner units. Two typical waves are presented. The first one, denoted "wave 1", is typical of a streaky structure (λ +

x ≈ 10λ + z ) representative of the near-wall cycle [START_REF] Jiménez | Near-wall turbulence[END_REF], for which the resolvent analysis is predictive (Abreu et al. 2020a,b). The second one, denoted "wave 2", has a similar phase speed, but is of larger size with a lower aspect ratio λ +

x ≈ 2λ + z . In the latter case, the turbulent structures are larger and extend further from the wall. For this type of waves, resolvent analysis shows significant discrepancies compared to the reference SPOD modes. Exact values of wavelengths are given in the figure captions. We have

α = 2π λ + x , β = 2π λ + z and ω = 2π λ + t
, where λ +

x and λ + z denote wavelengths in x and z directions, respectively, and λ + t is the time period, with all results expressed in inner units. Only antisymmetric waves (in u) with respect to the channel centerline are displayed. Such antisymmetric modes can be obtained by restricting the forcing in resolvent and stochastic models to be antisymmetric, and by taking the antisymmetric part of disturbances prior to computing SPOD. Similar analysis for symmetric waves shows similar results and will not be shown here for conciseness. Results are shown here for the u component, and corresponding results for v and w can be found in appendix B.

Considering wave 1 at Re τ = 180 in figures 2(a) and 2(c), ν-resolvent analysis is able to reproduce SPOD with a good accuracy. For this kind of waves, the improvement margin is quite low. In accordance with previous results in the literature, such as [START_REF] Morra | On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows[END_REF], incorporating the eddy viscosity in the resolvent analysis slightly improves further the agreement with SPOD. SLM with a noise based either on SPOD or resolvent modes have similar performances than ν t -resolvent to recover the first SPOD mode. In the Appendix, it is shown in figure 9 that for wave 1, SLM is slightly better to predict u than ν t -resolvent, and leads to a slight worsening for v and w. We can note a very good agreement of SLM with resolvent noise to predict u. The hierarchy of performances is less clear for mode 2, except that SLM with a SPOD noise improves significantly the prediction of u.

For wave 2 in figures 2(b) and 2(d), where ν-resolvent shows large discrepancies, the eddy viscosity improves well the behaviour. In that situation, SLM improves the agreement, especially for the noise based on resolvent analysis. In figure 10 in the Appendix, it can be shown that the improvement of v predictions is strong. This can be explained by the fact that the forcing shape equation (4.19) comes for the term coupling

+ x = 1124, λ + z = 102, λ + t = 100, Wave 2: λ + x = 2249, λ + z = 1124, λ + t = 200
. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), stochastic linearised model (SLM) with a noise defined by SPOD (blue) and resolvent (purple).

velocity components f x = -( σ ηα,β,ω ) y ∂U ∂y . The improvement is more pronounced for mode 2.

Comparison of the results in figures 2 and 3 shows that the performance of SLM at Re τ = 550 is quite similar to what is obtained for Re τ = 180. Exactly the same analysis can be done on wave 1. This gives us some confidence on the robustness of the modelling procedure. Wave 2 figure 3(b) is slightly different since it extends more in the center of the channel, from the lower Re τ results, as it represents larger structures extending to higher wall-normal locations, more relevant for larger Reynolds numbers. In this case, the disagreement between ν-resolvent analysis and SPOD is amplified. The ν t -resolvent finds a solution with a wider extend but located too close to the wall. SLM with SPOD seems to find a solution of the same nature than ν t -resolvent, but SLM with resolvent noise begins to have large values near the center of the channel as the SPOD. This supports the importance of modelling properly the noise statistics in the center of the channel.

Comparison of cross-spectral densities

For a more complete view, the CSD S α,β,ω (y, y ) = E ( ûα,β,ω (y) ûα,β,ω (y )) can be compared. This matrix is already computed to obtain SPOD (from the DNS data) and stochastic model (by the construction of an ensemble of realisations). As for the resolvent operator, assuming the system is excited by a Gaussian white noise, the cross spectral density can be predicted [START_REF] Cavalieri | Wave-packet models for jet dynamics and sound radiation[END_REF] by

S Resolvent α,β,ω = R α,β,ω E fα,β,ω f * α,β,ω R * α,β,ω = R α,β,ω R * α,β,ω .
(5.1)

Focusing on wave 1 at Re τ = 180, the cross spectral density matrix is compared in

λ + x = 1137, λ + z = 100, λ + t = 100, wave 2: λ + x = 3412, λ + z = 1706, λ + t = 200
. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), stochastic linearised model (SLM) with a noise defined by SPOD (blue) and resolvent (purple). figure 4. Only the streamwise component S uu α,β,ω (y, y ) = E (û α,β,ω (y)û α,β,ω (y )) is shown for sake of conciseness. It can be shown that ν-resolvent analysis predicts a too short coherence decay (seen by small values at off-diagonal terms in the CSD), with an elliptic pattern in the CSD, while the SPOD have rather a guitar plectrum shape. The prediction with ν t -resolvent is more accurate, but the non-ellipticity is too pronounced. Both SLM models improve significantly the prediction.

Collinearity metric

A comprehensive view of the modelling performances over all scales can be given by a collinearity criteria As a reference, figure 5 represents the collinearity between SPOD and eddy resolvent analysis for all wavelengths at λ + t = 200 and λ + t = 1000 for Re τ = 180. It shows that, in agreement with the analysis of turbulent pipe flow by Abreu et al. (2020b) (performed for ν-resolvent), the first ν t -resolvent mode predicts well the first SPOD modes for a wide range of wavelengths, especially for elongated modes (λ +

β model j,α,β,ω = Φ model j,α,β,ω , Φ SPOD j,α,β,ω Φ model j,α,β,ω Φ SPOD j,α,β,ω , (5.2
x > λ + z ), exemplified previously by wave 1. There is very little room for improvement in those regimes. The accuracy decreases for less elongated modes, of larger size and with a lower frequency (large λ + z and λ + t , low λ + x ), exemplified by wave 2. In all cases, the second ν t -resolvent mode leads poor predictions of the second SPOD mode. First SPOD Eigenvalue in figures 5(c) and 5(g) are shown to evaluate the presence of these structures in the flow in terms of kinetic energy. A similar behaviour is shown in appendix C for Re τ = 550 (figure 13). The worse agreement between SPOD and ν t -resolvent is directly associated with the decay of the ratio of the two first ν t -resolvent singular values (see figures 5(d) and 5(h)). Such cases require a finer description of the noise for accurate modelling. at Re τ = 180 and λ + t = 200. We can see that for elongated waves, where ν t -resolvent is highly performant, SLM has slightly lower agreement metrics, but there is a clear improvement to match SPOD at wavenumbers for which ν t -resolvent is less efficient.

In order to quantify the improvement (in terms of SPOD predictability) brought by the stochastic model compared with ν t -resolvent, we define the log-ratio of collinearity criteria γ SLM noise j,α,β,ω = log

β SLM noise j,α,β,ω β ν t -resolvent j,α,β,ω .
(5.

3)

The label "SLM noise" refers to the noise definition of the stochastic linearised model. It can be defined by SPOD (SLM SPOD) or resolvent (SLM Res.). Figures 7 and8 show γ SLM noise j,α,β,ω respectively at Re τ = 180 and Re τ = 550, for λ + t = 200. Other frequencies are displayed in appendix C (figures 14 to 17). In all these figures, red (blue) regions show parameters where SLM leads to better (worse) predictions compared to ν t -resolvent predictions.

As a first remark, for all Reynolds numbers and frequencies and for both noise definition, the trend described previously is confirmed, i.e. an improvement (γ SLM. noise j,α,β,ω > 0) at wavenumbers for which ν t -resolvent has moderate performances, and a slight performance loss for elongated waves. Moreover, the improvement is generally better for mode 2. Let us remark that for these elongated waves, the agreement is still good, and better than ν-resolvent analysis, as it has been shown by the profiles of wave 1. Considering the choice of the noise, SLM (SPOD) is a reference case, where the noise is defined from DNS data. Resolvent noise SLM (Res.) is a model-based technique whose goal is to approach the SPOD noise performance. In the latter case, the improvement follows the same trend than the one with SPOD noise, with a slightly lower performance (especially at Re τ = 550) than the one that can be expected since the modelling is now purely model-based. Resolvent modelling thus constitutes a robust and efficient way to define the noise. It clearly improves the performances of the resolvent analysis while remaining purely based on the model. Obviously, if SPOD modes are available, it is preferable to use them.

Discussion

The reason for the improvement of the present stochastic model compared to resolvent analysis is that SLM considers the effect of the turbulence on the waves by a stochastic transport, thus preserving conservation of quantities. Such effect is by essence non-Gaussian through the transport by the Wiener process. Moreover the stochastic diffusion, which may be interpreted as a turbulent viscosity, is constructed to dissipate the energy brought by the noise. This specific structure is a step to go beyond additive Gaussian noise that is implicit in resolvent analysis applied to turbulent flows.

In [START_REF] Symon | A tale of two airfoils: resolvent-based modelling of an oscillator versus an amplifier from an experimental mean[END_REF], an energy equation is written for individual resolvent modes, and energy transfers are discussed. The results in the cited work show that for frequencies and wavenumbers for which there is a large production, ν t -resolvent analysis is successful. Indeed, by its own nature the resolvent predicts well production, and the excess of energy that is not dissipated by molecular viscosity but drained by non-linear energy transfer is accurately modelled by eddy dissipation. However, for waves that gain energy from the others, the non-linear energy transfer cannot structurally be well represented by a purely diffusive mechanism. In SLM, the stochastic diffusion plays the role of eddy viscosity and is important to represent well waves for which high production occurs (in our case the elongated structures). In addition, the stochastic transport is a mechanism that can model energy transfers coming from other (unresolved) scales. Moreover, in the temporal formulation, these terms are designed such that energy is conserved. Obviously, this conservation is, by Parseval theorem, integrated over all frequencies. These energy transfer considerations explain the fact that SLM improves predictions for scales poorly modelled by the ν t -resolvent, thanks to the possibility of modelling gain of energy from other scales. We interpret the slight loss of performance for elongated scales as a room for improvement in the definition of the noise, so as to improve the properties of stochastic diffusion.

Noise definition in SLM can be viewed as a closure problem. The main improvement could be attributed to the fact that we have introduced data in the modelling of the background turbulence through SPOD. The fact that the improvement is still visible with the resolvent noise is a good indication that a true modelling improvement is obtained.

Conclusion

A stochastic procedure is proposed to model coherent structures, or waves, in turbulent flows. This procedure, which has been tested here for turbulent channel flow at two friction Reynolds numbers, is based on the resolution of an ensemble of linearised problems derived from a stochastic version of the Navier-Stokes equations. Assuming implicitly a triple decomposition, we consider a wave (low-amplitude solution of the model) evolving over a mean flow and perturbed by turbulent velocity fluctuations sustaining a stochastic transport. In the present framework, these fluctuations are modelled by a Wiener process, where spatial correlations are prescribed.

In order to define such stochastic properties of the background turbulence, two different noise definitions are proposed, based on SPOD and resolvent analysis: the former datadriven and the latter model-based which provide a first approximation of the overall background turbulence.

The background turbulence in the stochastic framework leads to a modified linearised operator, with an enhanced diffusion akin to an eddy viscosity, and a drift velocity that modifies the base flow. Moreover, the stochastic disturbances act as a forcing term. The background turbulence acts thus both injecting and dissipating energy of coherent waves.

In the present study, the drift velocity has been found to be mainly active in the buffer region. This is consistent with the study of [START_REF] Pinier | A stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF]. The stochastic diffusion is also largest at the buffer layer, where flow fluctuations peak.

When compared to predictions based on ν t -resolvent analysis, the present stochastic linearised model, with the two choices of noise modelling, shows improvement in the predictions of coherent structures in turbulent channel flow, for frequencies and wavenumbers for which the ν t -resolvent fails, with a slight performance reduction for waves for which the ν t -resolvent is highly accurate. This is quantified here by comparison with SPOD modes taken from direct numerical simulations. This behaviour is understood by the ability of SLM to directly model energy gains coming from other scales by non-linear transfer. The good performance shows the pertinence of the present linearised stochastic approach, which provides a consistent manner to model interactions between coherent structures and background turbulence. The linearisation implies that among all nonlinear interactions, the present strategy models the transport of the perturbation by the turbulent fluctuation.

Resolvent analysis employed under the simplistic white-noise forcing assumption leads to inaccuracies [START_REF] Nogueira | Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow[END_REF][START_REF] Morra | The colour of forcing statistics in resolvent analyses of turbulent channel flows[END_REF]). The present model shows a way for including interactions with the background turbulence, avoiding use of an ad hoc eddy viscosity, which, despite improved predictions in some cases [START_REF] Morra | On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows[END_REF], lacks generality. Of course, it is replaced by the definition of the statistics of a background turbulent field, with quantities possibly easier to educe or deduce than an eddy viscosity. Compared to strategies based on the stochastic forcing coloration, our method is based on velocity statistics, which are more direct to obtain and analyse than non-linear term statistics. We here provide a way to include interactions of coherent structure with background turbulence that is coloured in space but white in time, similar to what was previously used for estimation of turbulent channel flow [START_REF] Chevalier | State estimation in wall-bounded flow systems. Part 2. Turbulent flows[END_REF]. Such inclusion of colour explains the more robust predictions of the present model compared to resolvent analysis.

The noise definitions proposed in this paper are not the only possibilities, and a fine design of the associated statistics is an open research direction. SPOD and resolvent analysis are proposed here, but considerations based on turbulence modelling, or empirical Reynolds stresses (profiles often easily available) are promising potentialities. Statistics of the noise can be learned as well by formulating an inverse problem.

  y )φ j (y ) dy = λ j φ j (y) (3.2) leads to an orthonormal basis of SPOD modes (since S(y, y ) is Hermitian), with respect to the inner-product (a(y), b(y)) y = 1 -1 a(y)b(y) dy. The real eigenvalues λ j represent the relative contribution of a given mode to the CSD. In practice, SPOD modes are computed by splitting the signal into an ensemble of M possibly overlapping time windows of size N , whose snapshots u (n) i are spaced ∆t apart

1- 1

 1 σα,β (y, y ) ηα,β,ω (y )dy , with ηα,β,ω (y) = dη α,β,ω (y)/dt a standard centered Gaussian white noise profile.System (4.8) is linearised around the parallel (finite-time averaged) mean field ū = (U (y) 0 0) T . We define as well ūd = (U d (y) 0 0) T , with U d (y) = U (y) -1 2 ∂a xy ∂y . The drift contribution may be non-null. Indeed we have ∂a xy ∂y = ∂σ xx ∂y σ yx + σ xy ∂σ yy ∂y + σ xz ∂σ yz ∂y .

Figure 1 .

 1 Figure 1. Effect on the stochastic diffusion and the drift velocity at Re = 550.

Figure 2 .

 2 Figure 2. Profiles of power spectral density of |û| at Reτ = 180 for two waves. Wave 1:+ x = 1124, λ + z = 102, λ + t = 100, Wave 2: λ + x = 2249, λ + z = 1124, λ + t = 200. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), stochastic linearised model (SLM) with a noise defined by SPOD (blue) and resolvent (purple).

Figure 3 .

 3 Figure 3. Profiles of power spectral density of |û| at Reτ = 550 for three waves, Wave 1:λ + x = 1137, λ + z = 100, λ + t = 100, wave 2: λ + x = 3412, λ + z = 1706, λ + t = 200. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), stochastic linearised model (SLM) with a noise defined by SPOD (blue) and resolvent (purple).

  ) defined similarly as in Cavalieri et al. (2013); Lesshafft et al. (2019); Abreu et al. (2020b).

Figure 4 .

 4 Figure 4. Real part of S uu α,β,ω (y, y ) for wave 1 λ + x = 1137, λ + z = 100, λ + t = 100 at Reτ = 180.

Figure

Figure 5 .

 5 Figure 5. Value of β ν t -resolvent j,α,β,ω

  Mode 1, νt-resolvent.

Figure 6 .

 6 Figure 6. Value of β ν t -resolvent j,α,β,ω , β SLM (SPOD) j,α,β,ω

Figure 7 .

 7 Figure 7. Value of γ SLM noise j,α,β,ω , Reτ = 180, λ + t = 200.

  Mode 2, SLM (Res.).

Figure 8 .

 8 Figure 8. Value of γ SLM noise j,α,β,ω , Reτ = 550, λ + t = 200.

Figure 9 .

 9 Figure 9. Profiles of power spectral density at Reτ = 180 for wave 1: λ + x = 1124, λ + z = 102, λ + t = 100. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), Stochastic linearised model (SLM) with a noise defined by SPOD (blue) and resolvent (purple).

  Mode 2, | ŵ|.

Figure 10 .

 10 Figure 10. Profiles of power spectral density at Reτ = 180 for wave 2: λ + x = 3412, λ + z = 1706, λ + t = 200. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), Stochastic linearised model (SLM) with a noise defined by SPOD (blue) and resolvent (purple).

Figure 11 .

 11 Figure 11. Profiles of power spectral density at Reτ = 550 for wave 1: λ + x = 1137, λ + z = 100, λ + t = 100. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), Stochastic linearised model (SLM) with a noise defined by SPOD (blue) and resolvent (purple).

  Mode 2, | ŵ|.

Figure 12 .

 12 Figure 12. Profiles of power spectral density at Reτ = 550 for wave 2: λ + x = 3412, λ + z = 1706, λ + t = 200. SPOD (black), ν-resolvent modes (red), νt-resolvent (orange), Stochastic linearised model (SLM) with a noise defined by SPOD (blue), and resolvent (purple).

  νt-resolvent singular values for λ + t = 1000.

Figure 13 .

 13 Figure 13. Value of β ν t -resolvent j,α,β,ω

  Mode 2, SLM (SPOD).

  Mode 2, SLM (Res.).

Figure 14 .

 14 Figure 14. Value of γ SLM noise j,α,β,ω , Reτ = 180, λ + t = 500.

Figure 15 .

 15 Figure 15. Value of γ SLM noise j,α,β,ω , Reτ = 180, λ + t = 1000.

  Mode 2, SLM (SPOD).

  Mode 2, SLM (Res.).

Figure 16 .

 16 Figure 16. Value of γ SLM noise j,α,β,ω , Reτ = 550, λ + t = 500.

Figure 17 .

 17 Figure 17. Value of γ SLM noise j,α,β,ω , Reτ = 550, λ + t = 1000.

† Email address for correspondence: gilles.tissot@inria.fr

† It is abusively labelled as "resolved" as this component is of course specified by the user.

† More precisely, the model is bilinear with respect to (u , σdBt).

Declaration of interests

The authors report no conflict of interest.

Appendix A. The Stochastic Navier-Stokes equations

Mass conservation is obtained by applying the stochastic transport equation (4.6) to the density with θ = ρ. With incompressibility assumption, i.e. ρ = 1 (working with dimensionless variables), we obtain immediately

Momentum conservation is obtained with θ = ρu i with u i the i th component of u. This leads to

where F i (X t , t) is the i th component of the forces acting on the fluid particle. The term d t F i (X t , t) is the associated force impulse. Let us remark that

where ρ, u i is the quadratic variation between the two stochastic processes ρ and u i .

Rewriting in non-conservative form similarly than for the deterministic case and applying incompressibility condition leads to

Applied forces are the pressure and viscous stresses:

where p t is the finite-variation contribution of the pressure and dp t is a zero-mean stochastic process describing the pressure fluctuations due to the martingale part of the velocity component. Besides its physical meaning, this stochastic term is mandatory to balance the martingale part of the equations. Similarly, the viscous stresses induced by the fluid particle stochastic displacement are decomposed in a finite variation and martingale part. A rigorous justification of this decomposition based on conservation integrals is given in [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF].

Appendix B. Profiles of velocity fluctuations

In this section, amplitudes of the three velocity components are presented for waves 1 and 2 at both Reynolds numbers (figures 9 to 12). This complements the results presented in section 5 and shows that the stochastic modelling improves especially well u of wave 1 and v of wave 2.

Appendix C. Collinearity criteria

In this section complementary collinearity criteria are given. Figures 13 show the collinearity between ν t -resolvent modes and SPOD at Re τ = 550. It shows that all the trends remain similar to the ones observed at Re τ = 180.

Improvement of collinearity by the stochastic linearised model is shown for Re τ = 180 for λ + t = 500 (figure 14) and λ + t = 1000 (figure 15), and for Re τ = 550 for λ + t = 500 (figure 16) and λ + t = 1000 (figure 17). The behaviour is similar for all tested frequencies, and the improvement is more pronounced at Re τ = 180 than Re τ = 550.