Comptes Rendus Mécanique

. An attempt to improve the accuracy of resolvent-based predictions by including velocity correlations in the linear model is developed here. Closure assumptions for unresolved nonlinearities are thus pushed back to a higher order. Turbulent channel flow is considered as a test case: response and forcing modes obtained from singular value decomposition of the new resolvent model are compared to Spectral ProperOrthogonalDecomposition(SPOD)modesextractedfromaDirectNumericalSimulation(DNS)data-base. The performance of the approach is also measured against previous resolvent-based models. The new model does not yield significant global improvement, but does improve predictions in some regions. Further work on the method should target the linear modeling of the velocity-pressure gradient correlation tensor.


Introduction
Incompressible turbulent flows are very diverse and cover a broad range of applications.However, they are also notoriously difficult to model.In this work the choice was made to solve the Navier-Stokes equations (NSE) in frequency space instead of approaching time-dependent solutions directly.Here, a novel spectral reduced-order approach to model incompressible turbulent flows will be presented.This method is derived from resolvent modeling with an input from Reynolds-stress models.The resolvent approach proposes a linear model of a fully turbulent flow at a given frequency.Interpreting nonlinearities in the NSE as a forcing input into the linearized fluctuation equations, the most amplified response structures are taken as a prediction for energetic turbulent flow structures, for instance in [1,2].This approach is most often applied in frequency space and comes with all the advantages and drawbacks of a linear model, namely inability to account for all terms, simplicity of implementation, as well as computational efficiency.Reynolds-stress models originate from a theoretical attempt to improve the quality of Reynolds averaged Navier-Stokes models by adding second-order terms in the relevant equations as in [3, chapter 11].
Resolvent analysis has already demonstrated its ability to accurately predict large-scale structures in turbulent channel flow at a given frequency; [4] already exposed excellent agreement between Direct Numerical Simulation (DNS) data and a combination of twelve resolvent modes with weights derived from an optimisation problem computed at every wavenumber in a spectral framework.Said agreement persisted throughout the DNS energy spectrum and across the channel without any empirical modeling.Quantifying complexity in the same study clearly shows the main advantage of reduced-order modeling -indeed, the computations involved with the resolvent model are more than three orders of magnitude cheaper than the relevant DNS, and they nonetheless produce high quality estimations of the fluctuations.Despite these encouraging results, this approach provides no information on nonlinear terms in the original equations, which can also represent an engineering objective.
Indeed, other authors have given thought to this idea of a low-order modelling of turbulent flows that would be more cost-efficient than DNS -good enough for engineering purposes, and orders of magnitude less expensive.[5] consider a stratified channel flow and prove that on the limited range of scales considered, a resolvent model is able to reproduce an energy budget similar to that of a DNS.
As part of the effort to improve resolvent model performance for wall-bounded flows, [6] apply a decomposition of the velocity field in solenoidal and non-solenoidal parts in the Orr-Sommerfield-Squire (OSS) formulation.Exact Coherent States (ECS) for simple Couette and Poiseuille flows were used to validate the model's performance.An astute method to explicitly compute interaction coefficients as part of an optimisation problem led to a marked increase in performance.A single resolvent mode was demonstrated to accurately represent some ECS solutions.Finally, a characterization of nonlinear terms was performed using the DNS of a channel flow, but stopped at the conclusion that its spectral signature is non trivial.The most important drawback of this approach is that the OSS formulation makes the corresponding nonlinear terms much harder to interpret and model, thus limiting the possibility for further model improvement.
[7] computed Spectral Proper Orthogonal Decomposition (SPOD) and resolvent modes for a wide range of wavenumbers from a DNS database.That study avoids the use of eddy viscosity and presents large discrepancies between SPOD and resolvent modes.Possible explanations include different physical phenomena which dominate various flow regions, such as the lift-up effect which seems to produce favourable conditions for model and data agreement.[7, Figures 3 and  4] display the correlation between SPOD and resolvent modes which clearly illustrate the appeal of a new reduced-order method that would be more general than the current state of the art.[8] present good agreement between a DNS power spectral density (PSD) and a modeled one, using eddy viscosity and a resolvent formalism in a channel flow, as well as an engineered forcing that corresponds better to true flow statistics than white-noise.
More recently, [9] detailed how to obtain a "coloured" forcing term more representative of flow physics.This approach results in a tool with predictive value, though its generality is somewhat lessened by the use of an eddy viscosity model.
[10] present a resolvent-based estimation of turbulent channel flow using wall fluctuations as an input.This work also highlights the importance of "coloured" forcing in the NSE, as whitenoise forcing leads to large fluctuations in the logarithmic layer.
Another recent study by [11] stresses how the addition of eddy viscosity makes a significant difference in the resolvent model's ability to properly capture flow energy.Mainly concerned with the latter quantity, [11] outline the balance between production, dissipation, and nonlinear transfer of turbulent fluctuation kinetic energy.
Indeed, the ability of eddy viscosity to improve resolvent models and enable them to reproduce fine structures close to the wall present in DNS is desirable in channel flows, as stressed by [12].Although the addition of eddy viscosity makes the resolvent model more robust across wave speeds, the conclusion of [12] was that a scale dependant eddy viscosity should be considered.
A major drawback of resolvent methods listed above is their inability to predict nonlinear terms representative of the real flow.Indeed, the analysis of such a linear operator provides response modes of physical quantities and corresponding optimal forcing terms that represent flow nonlinearities.These forcing terms are rarely representative of actual physical flow quantities.In an effort to address this limitation, and in order to obtain the "coloured forcing" directly in the new approach, a usual resolvent model for the velocities was coupled with equations for velocity correlations such as to obtain a "second-order resolvent model" that would accurately account for both dominant order response and nonlinear terms in the NSE, called forcing terms.
This idea can to some degree be found in the stochastic structure stability theory of [13].Indeed, this work details how to obtain a linear model for a turbulent mean flow and the associated Reynolds stresses where both evolve in time.The authors acknowledge that for large Reynolds numbers the method becomes highly intractable, preferring instead a reduced order model which uses a few realisations of the fluctuations forced with random white noise to approximate non-linear terms in the baseflow equations as the system becomes large.
Another similar approach goes by the name of generalised quasilinear approximation.This formalism explicitly introduces a range of frequencies that are allowed to interact with one another.[14] explain how this is a generalisation of the quasilinear approach where a single frequency is considered independently of all others.Indeed, in order to close the model, some interactions between waves are cut off from the generalised model.This approach does not always fare well in describing the whole system, but almost always constitutes an improvement compared to the classical quasilinear approximation.
This paper is organised as follows: Section 2 presents the methodology of usual resolvent analysis as well as the additional equations specific to the new approach, Section 3 gives more details about parameter choices and the test dataset, Section 4 presents the two linear systems and links them together, Section 5 exhibits results across the wavenumber space and details a specific set of parameters.Finally, Section 6 summarises the study.

Methodology
This section aims at laying out necessary formalism and equations.Starting from the NSE, and a given baseflow, analytical tools to obtain insights into the coherent structures of the flow will be built.The first goal is to obtain a linear system u ′ = R f which transforms a given forcing f representing non-linear terms of the relevant equations into a fluctuation velocity field u ′ through a resolvent operator R .The standard approach of defining the linear resolvent operator, as widely applied in the recent literature, is presented in Section 2.1.The proposition for an extended linear resolvent operator is introduced in Section 2.2.Subsequently, Section 2.3 will explain how a linear system can be exploited to obtain its most amplified fluctuations and optimal forcing.

First-order equations
Writing the incompressible NSE for total velocity u = U +u ′ and pressure p = P +p ′ with ensemble averages u = U and p = P and the convention [∇U ] i j = ∂U i /∂x j gives with τ = u ′ u ′T , and the total derivative is defined as is nonlinear in the fluctuation quantities q u = [u ′T p ′ ] T , it may be interpreted as forcing as in [4].The linear model system for flow fluctuations is written Note that the first line of the linear system relates to the incompressibility condition, which is not forced.Resolvent analysis thus circumvents the closure problem by labelling unclosed terms as forcing.Conceptually, one goes from trying to solve for q u = [u ′T p] T to looking for coherent structures favoured by the linear operator L u .

Second-order equations
The key proposition of this paper is to model the τ tensor, in an attempt to increase the accuracy of the resolvent model.From equation (1) one obtains The right hand side of equation (3) will be considered as forcing terms, denoted f , their common trait being their nonlinearity with respect to u ′ , p ′ and τ .
Replacing the symmetric f and τ tensors with vectors f v and τ v , which each contain six independent tensor components, equation ( 3) can be cast as a linear system

Resolvent analysis
Let us define a generic extensor B , which extends a vector f with zeros Given a state vector q that includes the fluctuation velocities u ′ (for example q u = [u ′T p ′ ] T ), let a generic extractor (or observation matrix) H , be the operator that extracts the velocities from the entire state Assume a linear operator L that links the state vector q (which again includes the fluctuation velocities u ′ ) and the forcing vector f through an extensor B (as is the case in equation ( 2)) as follows L q = B f .
Using these general matrices B , H and L , which will be detailed later for specific systems, let the resolvent matrix be R = H L −1 B , so that one finally has u ′ = R f .The Singular Value Decomposition (SVD) of the R matrix, R = i σ (i ) ψ (i ) φ (i )H yields σ (i ) the gain for mode i of the R matrix, ordered in decreasing value.The normalised vectors φ (i ) (right singular vectors) thus computed are called in the following the forcing modes, and ψ (i ) (left singular vectors) will be referred to as the response modes of the resolvent system.The positive real-valued σ (i ) represents a gain between the norm of the forcing and that of the response.The first of these modes, φ (1) and ψ (1) , are optimal in the sense that they represent the forcing/response pair associated with the highest gain.This entire process is not new in the case of system (2) and was performed by [7,9], amongst others.

Flow configuration
The specific case of a fully developed channel flow between two infinite plates was studied, separated by a distance 2h.The presence of viscosity imposes u ′ = 0 at both walls.This naturally translates into τ = 0 at the wall.Even though there is physical motivation to impose additional boundary conditions for τ , the order of the system prevents the addition of further constraints.
For a wide range of parameters, forcing and response structures are concentrated near the channel walls, with very weak coupling between the top and the bottom halves of the channel.This weak coupling leads to poorly conditioned numerical systems when the full channel height is resolved.To address this problem, following [7], only the half-channel height was resolved, and symmetry conditions at the centre line imposed as follows.
With x streamwise, y wall-normal and z spanwise directions, one obtains U (x, y, z) = U (y)e x with U a function made strictly even across the channel.If u ′ z is an even (respectively odd) function, and thus symmetric (resp.anti-symmetric) across the channel centre line, then it follows from equation ( 9) that p ′ is an even (resp.odd) function.In turn, this gives odd (resp.even) u ′ y , and even (resp.odd) u ′ x .The parity of the six components of τ is constantu ′ x u ′ y and u ′ y u ′ z are always odd, u ′ x u ′ z is always even, regardless of the parity of u ′ z .For the sake of conciseness, the analysis was constrained to modes with even u ′ z , yet in the course of the study this choice had little influence on the results.Symmetry was enforced a posteriori, after DNS calculations and singular value decomposition, as in previous work by [7].In the following, only a half-channel is studied.

Fourier transforms
Consider Fourier transforms in space and time for all flow quantities q(x, y, z, t ) = q(y) exp (ik x x + ik z z − iωt ) (8) Note that all terms depending solely on averaged flow quantities disappear through Fourier transform at non-zero wavenumber.This is crucial because the forcing of equation ( 2) becomes f u = −∇ • τ and the link with system (4) is now apparent.
Writing δ j i the Kronecker delta, dU = dU /dy baseflow shear and z )/Re gives the channel flow case new directional formulations of equations ( 1) and (3),

Wall-normal discretisation
Discretisation of the full channel in the wall-normal direction was performed using n = 129 Chebyshev-Gauss-Lobatto nodes, with y = 0 being the centre of the channel, to allow for more points closer to the channel wall.The operator ∂ y is implemented to spectral order using Chebyshev polynomials following [15].A test case for convergence relative to n is detailed in Section 5.2.
In the following, the notation d will be used to denote operators discretised using this spacing.Consider an integration matrix W d so that y)dy in the chosen discretisation.Modes computed using the process detailed in Section 2.3 satisfy so they are orthonormal for an equal spacing discretisation only.This can be corrected as follows.
Using the Cholesky decomposition Then one can then take the Singular Value Decomposition (SVD) of the , and obtain the true modes of the system φ (i ) . These new modes are indeed orthonormal with respect to the chosen spacing,

DNS dataset
The analysis is based on spectral DNS data of a channel flow.The data was used to extract the time-averaged mean flow and to compute the leading SPOD modes following [16], against which the linear model results will be compared.This dataset has been validated in previous studies such as [7,9].Two Reynolds numbers are defined as where U bul k is the mean velocity, h is the channel half-width, ν is the kinematic viscosity, and u τ = max(dU )/Re is the friction velocity.Inner scalings are further defined as x + = x • u τ /ν and u + = u/u τ .Two DNS datasets are used, with (Re = 2 800, Re τ = 180) and (Re = 10 000, Re τ = 550).The time-wise Fourier transform necessary for SPOD computation, was performed using the Welch method with a second-order sinusoidal window of size N = 256, with 75% overlap.The Welch transform process also featured a force correction to counter the spurious influence of the windowing function as in [17].More information on the chosen DNS parameters can be found in Table 1, which presents the size of the box in streamwise and spanwise directions L x and L z , the grid steps in these directions in inner scaling ∆x + and ∆z + , as well as the smallest and largest steps in the normal direction ∆y + min and ∆y + max .

Final linear systems
In the previous two sections a variety of transformations and hypotheses that introduced much in the way of notation was detailed.For the rest of this study, the hat notation for all Fourier coefficients (see Subsection 3.2) will be removed, as well as the prime symbols on fluctuations (see Subsection 2.1), the d (see Subsection 3.3) and the v (see Subsection 2.2) superscript notations in the interest of clarity.In this section, the final form of the linear systems of interest will be detailed in Sections 4.1 and 4.2.Section 4.3 explains the procedure for the comparison between the linear model results and the leading SPOD modes of the DNS data.

First-order system
Equation 9 can be directly discretised to yield a 4n × 4n discrete system, where [dU ] The extensor B u and extractor H u operators therefore are non-square identity matrices, So far, this formulation is not novel.It has been done in [2,7,9,10].The process detailed in Sections 2.3 and 3.3 yields a resolvent matrix R 1 rst order for the so-called "first-order system", as well as associated response modes ψ (i ) 1 rst order and forcing modes φ (i ) 1 rst order .The currently most successful resolvent-based modelling of energetic coherent structures in channel flow uses a first-order formulation (Section 2.1) with an eddy viscosity model.For comparison, the Cess model employed in reference studies [9,12] is reproduced here using the friction Reynolds number as defined in equation (14).When applied, this model replaces the molecular diffusion term in equation ( 9), T diff,ν = ∆u/Re, by a more complete term T diff,ν t = ∇ • [(1 + ν t /ν)( ∇u + ∇u T )]/Re.In order to verify that the above model accurately represents the flow dynamics, a comparison between the mean flow computed using the eddy viscosity and the one obtained by DNS was performed.The computation of the mean flow from the eddy viscosity was done as in [18, section 2].The velocity profiles are compared in Figure 1.Differences between the two curves exist but remain small.This is the resolvent system used by [7,9], with an eddy viscosity representing at least part of the unknown nonlinear terms.Applying once more the process in Sections 2.3 and 3.3 yields a resolvent matrix R eddy viscosity , as well as associated response modes ψ (i )   eddy viscosity and forcing modes φ (i ) eddy viscosity .

Second-order system
The original proposition in the method is to include equation ( 10) inside a unified resolvent operator that includes velocity, pressure, and τ in the state vector.Recall that the motivation is to improve the capability of the resolvent operator to correctly predict flow physics both in terms of perturbations and their first moment.
Equation 10 translates to a 6n × 6n linear system with the vector notations τ and f (see Section 2.2) written as As before, f τ represents the forcing of the second order equations.Together, equations ( 15) and 18 form a general system formulation that is 10n × 10n in size, The equation above makes use of the coupling matrix A of size 4n × 6n defined as The extensor B and extractor H may be explicitly derived The process detailed in Sections 2.3 and 3.3 yields a new resolvent matrix R 2 nd order , with its associated response modes ψ (i ) 2 nd order and forcing modes φ (i ) . Note that the latter correspond to the forcing terms of equation ( 10) instead of those of equation 9.

Comparison of model versus data
The process detailed above was implemented in a Matlab code (version R2019b) using the svds (subset singular value decomposition) function library to obtain resolvent modes.
The most amplified response and forcing modes obtained by the resolvent model was compared with the most energetic SPOD modes obtained from DNS data.Said SPOD modes are computed using the method of snapshots detailed in [16].This method allows us to obtain SPOD velocity modes, denoted ξ (i ) .The same approach was used to compute SPOD nonlinear modes, obtained directly from the nonlinear term of equation ( 1) in the DNS.SPOD modes represent the orthonormal basis that captures the most of the flow energy for any given rank as established in [2].They are obtained as eigenvectors of a Cross-Spectral Density (CSD) matrix computed at a given frequency.SPOD may be performed for both velocities and nonlinear terms of the NSE.While the gain-based ranking of the modes gives an indication of each mode's prevalence in the flow, there is no causal association between SPOD nonlinear and forcing modes obtained in this manner.In other words, ξ (i ) is not related to ζ (i ) , and both modes are computed independently.
The matrix system detailed in equation ( 19) is exact but the flow response q can only be determined if the forcing f u is known.Because of the large gain σ (1) and large gain separation σ (1) /σ (2) typically observed in a resolvent operator, it is expected that the first left-singular vector ψ (1) approximates ξ (1) , the actual flow fluctuations obtained from SPOD.This has been observed in previous studies such as [7,9], and the coherent structures in a flow are of engineering and scientific interest by themselves, especially if they are provided by a numerically cheap linear analysis process.
However, as pointed out for instance by [19, section III.c] and by [20], there is no reason why the nonlinear terms in the actual flow ζ (i ) should be approximated by the right-singular vectors φ (i ) 1 rst order of the resolvent matrix R 1 rst order .Indeed ζ (1) represent the best low-rank approximation of flow nonlinear terms, whereas φ (1)   1 rst order only represents the optimal input to obtain ψ (1)   1 rst order with respect to R 1 r st or d er .
Thus, there is no real physical insight to explain why ζ (1) and φ (1)   1 rst order would have a similar structure, yet it is accurate to say that only the part that is strongly amplified by the linear operator is expected to have an important influence on flow behaviour.By contrast, ∇ • R τ φ (1)   2 nd order includes additional physics that might yield a better approximation quality of ζ (1) .

Parameter space exploration
For a first comparison with results from [7,9], the agreement between SPOD data and resolvent predictions was evaluated by varying λ + = [λ + x λ + z ] T .The + superscript denotes inner scaling as defined in Section 3.4.This agreement is quantified in terms of the projection coefficient β = |ψ (1) W ξ (1)H |, i.e. between first resolvent mode and first SPOD velocity mode.As stated in Section 4.3, the dominant velocity mode obtained through SPOD from DNS data is compared with a left-singular vector of the resolvent operator.
Only the Re τ = 180 (Re = 2 800) case will be discussed here, but the Re τ = 550 (Re = 10 000) case was also explored with similar conclusions.Isocontours of β for the range of λ + available at the fixed frequency f + = 10 −2 are shown in Figure 2a for the new model and Figure 2b for the eddy viscosity enhanced first model.This graph indicates no clear incentive to use the "second-order" resolvent model over the first-order eddy viscosity resolvent model throughout the λ + plane at this frequency.Even though the new model performs slightly better over a restricted area in the low λ + x , high λ + z regime, its performance drops elsewhere.The value f + = 10 −2 was chosen as representative of the near-wall cycle, for consistency with [7].
Results obtained with different frequencies (10 −4 , 10 −3 , 10 −1 , 1) show a similar pattern.The isocontours are also shown at f + = 10 −4 in Figures 2c and 2d, where one can see that the eddy viscosity outperforms the second-degree resolvent model throughout the domain overall.

Near wall structures
In the following, a case study for (λ ) is presented.This choice of parameters corresponds to the "near-wall" case discussed in [9, paragraph 3.1].Indeed, these wavenumbers represent a particularly energetic point in frequency space for y + ≈ 15, so relatively close to the wall.
To verify that results are independent of the number of points n, an additional set of calculations at n = 257 was performed for the near wall regime.This calculation with roughly double the y resolution as the one presented below differed by less than 2% from the following results in all aspects.Therefore, it is considered proven in the rest of this study that the results are converged with respect to n.
Figure 3 shows the absolute value of the first SPOD velocity and resolvent response modes for a variety of models along the three flow directions.It features the "first-order model" from Section 4.1 with and without eddy viscosity, as well as the "second-order model" detailed in Section 4.2.Likewise, Figure 4 presents the first SPOD nonlinear and resolvent forcing modes for all directions.As detailed in Sections 2 and 4.3, the modes compared here are the first modes with respect to gain obtained from eigenvalue decomposition of a CSD matrix, and from the SVD of a resolvent matrix for a specific model.There is a noteworthy exception to this, "second-order forcing terms" computed as normalised ∇ • R τ φ (1)   2 nd order .Thus, all modes presented below verify x H W x = 1.As discussed in Subsection 3.1, only half a channel was considered supposing the velocity along z to be even for all modes.
All in all, Figure 3 yields overall positive results.In the streamwise direction (Figure 3a), where the fluctuations are most intense, the "second-order" resolvent model succeeds in capturing most of the SPOD mode.The model captures very well close wall behaviour in this direction with no need for empirical eddy viscosity or wall functions.However, the new resolvent model does tend to underestimate the fluctuation amplitudes in the other two directions (Figures 3b and 3c), performing rather poorly in the spanwise direction, where it smoothens out the flow oscillations to an extreme degree.As in [9], the introduction of an eddy viscosity model greatly improves the performance of the resolvent model (consider the difference between the eddy viscosity and the first-order curves).
Figure 4 shows that the nonlinear terms of equation ( 1) are better represented overall along the three directions using the second-order model than the first-order one.This result was to be expected as there is no reason to expect that φ (1)   1 rst order ≈ ζ (1) (see Section 4.3).Notice though that the forcing components in the wall-normal and spanwise directions are strongly underpredicted by the new resolvent model.The new model allows again for an accurate restitution of actual flow statistics close to the wall in the streamwise direction.The capability to predict, to some degree, the properties of the nonlinear terms may be useful to construct estimators, as proposed in [10].

Large scale structures
In the following, a case study for (λ x = 4.19, λ z = 1.26, c + = 16) is explored.As before, this choice of parameters and the name of this test case correspond to that of [9, paragraph 3.1].The denomination "large scale" is inferred from the energetically significant dynamics associated with these wavenumbers, which exhibit structures that span the entire channel.The scaling is identical to the previous section, and only the absolute value of the modes is presented.Figure 5 presents the first SPOD velocity and resolvent response modes, whereas Figure 6 presents the first SPOD nonlinear and resolvent forcing modes.
Overall, this case study yields more mixed results.In two of the three directions, the model performs worse than the basic first-order resolvent model or the eddy viscosity system.The new model consistently underestimates real statistics for low y + , and does not improve predictions compared to the simpler first-order model, even in the dominant direction.Worse, the secondorder model smoothens out variations of velocity amplitude in the spanwise direction.Here again, the addition of an eddy viscosity yields significant improvements in the predictions from the first-order resolvent model.
Interestingly, the model overestimates the forcing in the x direction in Figure 6 and underestimates the rest.In comparison, an eddy viscosity model yields somewhat relevant predictions in the spanwise direction but otherwise fails to capture the dynamics of the nonlinear terms in the flow.Again, this is to be expected, as optimal forcing of the resolvent matrix cannot be expected to correspond to actual flow statistics.Comparison between the first SPOD mode of the DNS data ξ (1) and predictions from various resolvent models ψ (1) .Velocity components are shown for (Re τ , λ x , λ z , c + ) = (180, 4.19, 1.26, 16).   (and prediction from second-order resolvent analysis ∇ • R τ φ (1) .Optimal forcing obtained from first-order resolvents are also shown for reference.Nonlinear and forcing components are shown for (Re τ , λ x , λ z , c + ) = (180, 4.19, 1.26, 16).

Conclusion
The resolvent operator, which relates forcing input of a linear system to the associated response, has been formulated for the Navier-Stokes equations such that the forcing input includes thirdorder turbulent correlations, which drive a velocity response indirectly via the correlation tensor τ .In contrast, classical resolvent analysis, classified as "first-order" in this paper, defines forcing input as a substitute for second-order correlations, which drive the velocity response directly.This study tested whether the shift of forcing to higher order improves the linear response predictions in the case of turbulent plane channel flow.The relevant quantities to compare are the leading SPOD modes extracted from turbulent DNS data, which represent energetic coherent turbulence structures, and the leading SVD modes from the linear resolvent model.In the present case of a flow with streamwise and spanwise invariant statistics, such modes are obtained for a range of wavenumbers for a logarithmic range of frequencies.Reasoning was done on a halfchannel and presented results for an even u z .
Classical first-order resolvent analysis generally yields meaningful predictions of the leading SPOD modes, but with much room for improvement in quantitative accuracy.The limiting assumption of this approach is that the nonlinear terms in the fluctuation equations can be replaced by white noise.The currently best available strategy is to use this assumption in combination with the inclusion of an eddy viscosity in the linear operator.Much recent discussion in the community revolves around the question of how to improve the predictions further by replacing the white noise with "coloured" noise input.The present approach instead aims to produce the "colour" of ∇ • τ by computing the τ in a spectral setting as a response to white noise forcing at a higher order.
The success of this strategy, as far as it has been possible to establish, is not as clear as might have been hoped.Streamwise velocity fluctuation amplitudes are very well predicted in a "nearwall" setting in the sense of [7,9], and the streamwise component of ∇ • τ obtained from the model also compares well to the corresponding SPOD mode from the DNS.On the other hand, wall-normal as spanwise velocity and ∇• τ amplitudes are strongly underpredicted.In the "largescale" setting in the sense of [7,9], our second-order model clearly gives less accurate results than the standard first-order resolvent, especially when eddy viscosity is included in the latter.This observation stands even if an imperfect model is used for the eddy viscosity, namely one that produces a baseflow that differs slightly from the DNS one.When the model performance is measured in terms of projection with SPOD modes throughout the parameter space, the secondorder model cannot compete with the first-order model with eddy viscosity.
The roadblock to accurate second-order resolvent modelling is a typical closure problem, probably caused by the velocity-pressure-gradient correlation detailed in Section 2.1.Efforts to cast this tensor in the form of a linear function of u, p or τ so that it could be incorporated into the resolvent operator, have not been fruitful.Yet the present study represents a first step in a new direction for improved resolvent analysis applied to turbulent flow, from where new paths are hoped to emerge.

Figure 1 .
Figure 1.Mean velocity profiles from DNS and computed as a steady state with the Cess eddy viscosity model (equation (17)).

Figure 6 .
Figure 6.Comparison between the first SPOD mode of ∇ • τ in the DNS data ζ (1) and

Table 1 .
Characteristic parameters of the DNS database.