# On the use of eddy viscosity in resolvent analysis of turbulent channel flow

Sean Symon<EMAIL_ADDRESS>Aerodynamics and Flight Mechanics, Faculty
of Engineering and Physical Sciences, University of Southampton, SO17 1BJ,
United Kingdom  Anagha Madhusudanan  Department of Applied Mathematics and
Theoretical Physics, Centre for Mathematical Sciences, University of
Cambridge, Cambridge CB3 0WA, United Kingdom  Simon J. Illingworth  Department
of Mechanical Engineering, University of Melbourne, Parkville, VIC 3010,
Australia  Ivan Marusic  Department of Mechanical Engineering, University of
Melbourne, Parkville, VIC 3010, Australia

###### Abstract

The predictions of resolvent analysis for turbulent channel flow are evaluated
for a friction Reynolds number of $Re_{\tau}=550$. In addition to the standard
resolvent operator where only the kinematic viscosity appears, a resolvent
operator augmented with the Cess eddy viscosity profile is considered. Adding
eddy viscosity significantly alters the low-rank behavior of the resolvent.
Regardless of the wave speed selected, the eddy resolvent is low-rank for
spanwise wavelengths of $\lambda^{+}_{y}=80$ and $\lambda_{y}/h=3.5$ in
comparison to the standard resolvent whose low-rank behavior depends on the
wave speed. The leading eddy resolvent modes are shown to have higher
projections onto the leading mode from spectral proper orthogonal
decomposition in comparison to standard resolvent modes. Neither analysis,
however, reliably predicts the most energetic wave speed. The standard
resolvent tends to overestimate it while the eddy resolvent underestimates it.
For scales where the most energetic wave speed is underestimated, the eddy
resolvent modes are energetic too close to the wall. The eddy resolvent does,
however, correctly identify the most energetic wave speed and mode shapes for
structures that are associated with the near-wall cycle or that are most
energetic at $z/h=\pm 0.5$. It is argued that these types of structures are
likely to be correctly predicted for any friction Reynolds number due to the
inner and outer scaling of the Cess eddy viscosity profile. Finally, it is
shown that the accuracy of eddy predictions relies on striking the right
balance between positive and negative energy transfers. Even though the eddy
viscosity primarily adds dissipation, its wall-normal gradient injects energy
in the near-wall region, resulting in mode shapes that are “attached” to the
wall. For some scales, however, the predicted positive energy transfer is too
strong thus biasing structures towards the wall. The ability of the Cess eddy
viscosity profile to model both positive and negative energy transfers
suggests that it could be optimized for individual scales to provide better
low-order models of turbulent channel flow.

††preprint: APS/123-QED

##  I Introduction

The Navier-Stokes equations linearized around the mean (time-averaged) flow
have been used to identify coherent structures in a variety of flows. In the
resolvent-based approach of Ref. McKeon and Sharma (2010), for example, the
linearized Navier-Stokes equations are analyzed from an input-output
perspective. The input is made up of nonlinear perturbation terms that are
treated as an intrinsic forcing to the linear resolvent operator and the
output is the perturbation velocity field. In the context of wall-bounded
flows, resolvent analysis has been exploited for a variety of applications
from finding reduced-order models of exact coherent states Sharma  _et al._
(2016) to predicting statistics of high-Reynolds-number turbulence Moarref
_et al._ (2013); Skouloudis and Hwang (2021). The ability of resolvent
analysis to identify prominent linear mechanisms, moreover, has made it an
attractive alternative to direct numerical simulation (DNS) for designing flow
control strategies Luhar  _et al._ (2014); Nakashima  _et al._ (2017); Toedtli
_et al._ (2019).

From DNS and experimental studies of wall-bounded flows, it is known that the
production of turbulent kinetic energy is mainly driven by the exchange of
energy from the mean flow to the fluctuations. Resolvent analysis models well
this exchange from the mean to the fluctuations, where the mean is assumed to
be known a priori. Therefore, the term in the energy budget that the resolvent
analysis captures most successfully is production Symon  _et al._ (2021); Jin
_et al._ (2021). The resolvent operator also tends to be low-rank for energy-
producing scales Moarref  _et al._ (2013). The term that resolvent analysis
captures the least accurately for an arbitrary scale is nonlinear transfer
between scales. One way to model this term is to add an eddy viscosity to the
resolvent operator Reynolds and Hussain (1972); Pujals  _et al._ (2009); Hwang
and Cossu (2010); Morra  _et al._ (2019); Kuhn  _et al._ (2021); Morra  _et
al._ (2021); Pickering  _et al._ (2021). The eddy viscosity provides
additional dissipation that removes energy from all scales. As such, it
attempts to model this nonlinear transfer of energy from the large scales to
the small scales.

In terms of structures, many studies have observed good agreement between the
structures predicted by resolvent analysis, and those found in DNS. The
leading resolvent mode, which is computed as the leading left singular vector
of the resolvent operator, is the dominant structure predicted by resolvent
analysis. The dominant mode from DNS is computed as the leading mode from
spectral proper orthogonal decomposition (SPOD) Lumley (1967); Picard and
Delville (2000). In fact, if the nonlinear forcing is white in space and time,
then resolvent and SPOD modes are theoretically equivalent Towne  _et al._
(2018). The nonlinear forcing, however, is not white in space and time Morra
_et al._ (2021); Nogueira  _et al._ (2021); Amaral  _et al._ (2021). As such,
the role of eddy viscosity is to model the effect of the nonlinear forcing
such that white in space and time forcing is sufficient to predict the correct
structures. Reference Morra  _et al._ (2019) showed that adding eddy viscosity
improved predictions of coherent motions in turbulent channel flow at a
friction Reynolds number of $Re_{\tau}=1007$. However, they only considered
the streamwise velocity component of two of the energetic scales in the flow.
Reference Symon  _et al._ (2020) compared resolvent predictions with and
without eddy viscosity for turbulent channel flow at $Re_{\tau}=2003$. The
study was limited, however, to temporal snapshots of the largest scales. The
most comprehensive comparison between resolvent and SPOD modes has been
presented by Ref. Abreu  _et al._ (2020), who computed the projection of the
leading resolvent mode onto the leading SPOD mode for turbulent pipe flow at
low Reynolds numbers. The authors noted good agreement, i.e. high projections
between the resolvent and SPOD modes, for scales where the lift-up mechanism
Ellingsen and Palm (1975); Landahl (1980) was active. They did not, however,
consider the improvements that could be achieved by adding an eddy viscosity
model as done by Ref. Tissot  _et al._ (2021).

A rigorous comparison of resolvent and SPOD modes across all scales and wave
speeds, therefore, is needed for resolvent analysis with and without eddy
viscosity. As such, one of the principal objectives of this article is to
quantify the projections of the leading resolvent mode computed with and
without eddy viscosity onto the leading SPOD mode. Previous work Symon  _et
al._ (2020) suggests that analyzing the accuracy with which the resolvent
model predicts the wave speed at which a structure convects is crucial when
comparing the model predictions with DNS. This motivates the second objective
of this work which is to perform a detailed analysis of the most amplified
wave speeds identified by resolvent analysis with and without eddy viscosity
and how they compare to the most energetic wave speed computed from DNS. The
most energetic wave speed of a structure is also indicative of the wall-normal
location where most of its energy is concentrated and, therefore, of the shape
of the mode in the wall-normal direction. The final objective of this article
is to provide insight into the types of nonlinear interactions that the eddy
viscosity is attempting to model. As mentioned earlier, the eddy viscosity
primarily adds extra dissipation to the energy balance, but its wall-normal
gradient has been shown to add energy overall to a scale Symon  _et al._
(2021). It is worth examining the wall-normal profiles of the transfers
introduced by eddy viscosity to determine if they model the positive energy
transfer in the near-wall region that has been observed in DNS for large
scales Lee and Moser (2019); Kawata and Tsukahara (2021); Doohan  _et al._
(2021); Hernández  _et al._ (2021).

The rest of the paper is organised as follows. Section II presents the
governing equations for channel flow and provides a brief overview of
resolvent analysis, the Cess eddy viscosity profile, and SPOD. The details of
the DNS for $Re_{\tau}=550$ are provided in Sec. III. The predictions of
resolvent and eddy analysis are evaluated in Sec. IV using scalar quantities.
These include the ratio of energy contained in the first pair of singular
values compared to the total energy, the projections of the leading resolvent
modes onto the leading SPOD modes, and comparing the most amplified wave speed
from resolvent analysis to the most energetic wave speed in DNS. Section V
provides a more detailed comparison of the mode shapes for select scales. The
extent to which a constant eddy viscosity profile, i.e. one that does not
depend on space, can successfully predict structures is also investigated in
Sec. V. The energy transfer processes introduced by eddy viscosity are
discussed in Sec. VI. The role of the eddy viscosity gradient, in particular,
is examined and artificially manipulated to understand its influence on the
mode shapes. Finally, conclusions and implications for optimizing a scale-
dependent eddy viscosity are suggested in Sec. VII.

##  II Methodology

Section II.1 describes the governing equations for plane Poiseuille flow and
their non-dimensionalization. A brief overview of resolvent analysis is
provided in Sec. II.2. In Sec. II.3, a modified resolvent operator, which
includes the Cess eddy viscosity profile, is formulated. In order to assess
the predictive capability of resolvent analysis with and without eddy
viscosity, the leading modes are compared to SPOD modes, which are computed
using the procedure summarized in Sec. II.4.

###  II.1 Plane Poiseuille flow equations

The non-dimensional Navier-Stokes equations for statistically steady,
turbulent plane Poiseuille flow are

$\frac{\partial\bm{u}}{\partial
t}+\bm{u}\cdot\bm{\nabla}\bm{u}=-\bm{\nabla}p+\frac{1}{Re_{\tau}}\bm{\nabla}^{2}\bm{u},$
(1a) $\bm{\nabla}\cdot\bm{u}=0,$ (1b)

where $\bm{u}(\bm{x},t)=[u,v,w]^{T}$ is the velocity in the $x$ (streamwise),
$y$ (spanwise), and $z$ (wall-normal) directions, $p(\bm{x},t)$ is the
pressure, and $\bm{\nabla}=[\partial/\partial x,\partial/\partial
y,\partial/\partial z]^{T}$. The friction Reynolds number
$Re_{\tau}=u_{\tau}h/\nu$ is defined in terms of the friction velocity
$u_{\tau}$, channel half height $h$, and kinematic viscosity $\nu$. Periodic
boundary conditions are applied in the streamwise and spanwise directions and
no-slip boundary conditions are imposed at the walls. The velocities are non-
dimensionalized by $u_{\tau}$, the spatial variables by $h$, and the pressure
by $\rho u_{\tau}^{2}$ where $\rho$ is the density of the fluid. A ‘$+$’
superscript denotes spatial variables that have been normalized by the viscous
length scale $\nu/u_{\tau}$.

###  II.2 Resolvent analysis

Equation (1) is Reynolds-decomposed leading to the following equations for the
fluctuations:

$\frac{\partial\bm{u}^{\prime}}{\partial
t}+\bm{U}\cdot\bm{\nabla}\bm{u}^{\prime}+\bm{u}^{\prime}\cdot\bm{\nabla}\bm{U}+\bm{\nabla}p^{\prime}-\frac{1}{Re_{\tau}}\bm{\nabla}^{2}\bm{u}^{\prime}=-\bm{u}^{\prime}\cdot\bm{\nabla}\bm{u}^{\prime}+\overline{\bm{u}^{\prime}\cdot\bm{\nabla}\bm{u}^{\prime}}=\bm{f}^{\prime},$
(2a) $\bm{\nabla}\cdot\bm{u}^{\prime}=0,$ (2b)

where $(\overline{\cdot})$ and $(\cdot)^{\prime}$ denote a time-average and
fluctuation, respectively. The mean velocity profile $\bm{U}=[U(z),0,0]^{T}$
is assumed to be known a priori from DNS. Equation (2) is written such that
all linear terms appear on the left-hand side. The nonlinear terms on the
right-hand side are lumped together as a forcing $\bm{f}^{\prime}$. Equation
(2) is Laplace-transformed in time and Fourier-transformed in the homogeneous
directions $x$ and $y$

$\hat{\bm{u}}(k_{x},k_{y},s)=\frac{1}{(2\pi)^{3}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\bm{u}^{\prime}(x,y,z,t)e^{st-
ik_{x}x-ik_{y}y}dxdydt.$ (3)

Upon integration of Eq. (3), we set $s=i\omega$ to consider the frequency
response $\hat{\bm{u}}(\bm{k})$ where $(\hat{\cdot})$ denotes the Fourier-
transformed coefficient and the wavenumber triplet
$\bm{k}=(k_{x},k_{y},\omega)$ consists of streamwise wavenumber $k_{x}$,
spanwise wavenumber $k_{y}$, and temporal frequency $\omega$. The equivalent
wavelengths in the streamwise and spanwise directions are
$\lambda_{x}=2\pi/k_{x}$ and $\lambda_{y}=2\pi/k_{y}$. The wavenumbers are
non-dimensionalized by $(1/h)$ and the wavelengths by $h$.

Equation (3) is substituted into Eq. (2) and rearranged into state-space form
(Jovanović and Bamieh, 2005)

$i\omega\hat{\bm{q}}(\bm{k})=\bm{A}(k_{x},k_{y})\hat{\bm{q}}(\bm{k})+\bm{B}(k_{x},k_{y})\hat{\bm{f}}(\bm{k}),$
(4a) $\hat{\bm{u}}(\bm{k})=\bm{C}(k_{x},k_{y})\hat{\bm{q}}(\bm{k}),$ (4b)

where the state $\hat{\bm{q}}$ consists of the wall-normal velocity $\hat{w}$
and wall-normal vorticity $\hat{\eta}=ik_{y}\hat{u}-ik_{x}\hat{v}$. The
matrices $\bm{A}$, $\bm{B}$, and $\bm{C}$ are the discretized forms of the
linearized Navier-Stokes operator, the forcing operator and the output
operator, respectively, and are defined in the Appendix. It is worth noting
that $\bm{A}$, $\bm{B}$, and $\bm{C}$ are independent of $\omega$ but are
functions of the wavenumber pair $(k_{x},k_{y})$ under consideration. For the
sake of brevity, this dependence is omitted for the rest of the paper.

Equation (4) is recast into input-output form

$\hat{\bm{u}}(\bm{k})=\bm{C}(i\omega\bm{I}-\bm{A})^{-1}\hat{\bm{f}}(\bm{k})=\mathcal{H}(\bm{k})\hat{\bm{f}}(\bm{k}),$
(5)

where $\mathcal{H}(\bm{k})$ is a linear operator called the resolvent that
relates the input forcing $\hat{\bm{f}}(\bm{k})$ to the output velocity
$\hat{\bm{u}}(\bm{k})$. Even if $\hat{\bm{f}}(\bm{k})$ is unknown, the
resolvent operator can be characterized by the singular value decomposition

$\mathcal{H}(\bm{k})=\hat{\bm{\Psi}}(\bm{k})\bm{\Sigma}(\bm{k})\hat{\bm{\Phi}}^{*}(\bm{k}),$
(6)

where
$\hat{\bm{\Psi}}(\bm{k})=[\hat{\bm{\psi}}_{1}(\bm{k}),\hat{\bm{\psi}}_{2}(\bm{k}),\cdots,\hat{\bm{\psi}}_{p}(\bm{k})]$
are the resolvent modes, which form an orthogonal basis for velocity, and
$\hat{\bm{\Phi}}(\bm{k})=[\hat{\bm{\phi}}_{1}(\bm{k}),\hat{\bm{\phi}}_{2}(\bm{k}),\cdots,\hat{\bm{\phi}}_{p}(\bm{k})]$
are the resolvent forcing modes which form an orthogonal basis for the
nonlinear forcing. $\bm{\Sigma}(\bm{k})$ is a diagonal matrix that ranks the
$p$th structure by its gain $\sigma_{p}(\bm{k})$ using an inner product that
is proportional to its kinetic energy, i.e.
$\left<\hat{\bm{\psi}},\hat{\bm{\psi}}\right>=\int_{-h}^{h}\hat{\bm{\psi}}^{*}\cdot\hat{\bm{\psi}}dz$.
The structure $\hat{\bm{\psi}}_{1}(\bm{k})$ is, therefore, referred to as the
optimal or leading resolvent mode and is the most amplified response of the
linear dynamics contained in the resolvent.

The true velocity field from experiments or DNS can be expressed as a weighted
sum of resolvent modes

$\hat{\bm{u}}(\bm{k})=\sum_{p=1}^{N}\hat{\bm{\psi}}_{p}(\bm{k})\sigma_{p}(\bm{k})\chi_{p}(\bm{k}),$
(7)

where $\chi_{p}(\bm{k})$ is the projection of $\hat{\bm{\phi}}_{p}(\bm{k})$
onto $\hat{\bm{f}}(\bm{k})$, i.e.

$\chi_{p}(\bm{k})=\left<\hat{\bm{f}}(\bm{k}),\hat{\bm{\phi}}_{p}(\bm{k})\right>.$
(8)

It can be noted from Eqs. (7) and (8) that if $\hat{\bm{f}}(\bm{k})$ is white
noise or projects equally onto the resolvent forcing modes, then the
contribution of a resolvent mode $\hat{\bm{\psi}}_{p}$ in reconstructing the
velocity field according to Eq. (7) is solely dependent on the associated
singular value $\sigma_{p}$. Moreover, if $\sigma_{1}\gg\sigma_{2}$, then it
is often argued that the velocity response can be well-approximated by the
leading resolvent response mode $\hat{\bm{\psi}}_{1}$ alone.

###  II.3 Cess eddy viscosity model

It has been shown, however, that the nonlinear forcing may have little to no
overlap with the leading forcing mode Morra  _et al._ (2021); Symon  _et al._
(2021); Barthel  _et al._ (2021) resulting in $\chi_{1}\ll\chi_{p\neq 1}$.
This nonalignment stems from the fact that the resolvent does not model well
the inter-scale nonlinear transfer. To address this shortcoming, an eddy
viscosity can be added to the linearized Navier-Stokes equations after
performing a triple decomposition of the total velocity field $\tilde{\bm{u}}$
into a mean component $\bm{U}$, coherent motions $\bm{u}$, and incoherent
fluctuations $\bm{u}^{\prime}$ Reynolds and Hussain (1972). A new set of
equations govern the coherent velocity and pressure:

$\frac{\partial\bm{u}}{\partial
t}+\bm{U}\cdot\bm{\nabla}\bm{u}+\bm{u}\cdot\bm{\nabla}\bm{U}+\bm{\nabla}p+\nabla\cdot\left[\nu_{T}(\bm{\nabla}\bm{u}+\bm{\nabla}\bm{u}^{T})\right]=\bm{d},$
(9)

where $\nu_{T}(z)$ is the total effective viscosity and
$\bm{d}=-\bm{u}\cdot\bm{\nabla u}+\overline{\bm{u}\cdot\bm{\nabla}\bm{u}}$ is
the disturbance term. Similar to Refs. Reynolds and Tiederman (1967) and Hwang
and Cossu (2010), the Cess Cess (1958) eddy viscosity profile

$\nu_{T}(z)=\frac{\nu}{2}\left(1+\left[\frac{\kappa}{3}(1-z^{2})(1+2z^{2})(1-\text{exp}(|z-1|Re_{\tau}/A)\right]^{2}\right)^{1/2}+\frac{\nu}{2},$
(10)

is employed in this study. The constants $\kappa=0.426$ and $A=25.4$ are
chosen based on a least-squares fit to experimental mean velocity profiles at
$Re_{\tau}=2000$ del Álamo and Jiménez (2006). Even though the Reynolds number
in this study is lower than the Reynolds number for which the fit was
performed, it has been verified that the results are not sensitive to the
values of these constants.

As done in Sec. II.2, Eq. (9) is Fourier-transformed in time and the
homogeneous directions to obtain an input-output relationship between the
velocity and disturbance fields

$\hat{\bm{u}}(\bm{k})=\mathcal{H}^{e}(\bm{k})\hat{\bm{d}}(\bm{k}),$ (11)

where $\mathcal{H}^{e}(\bm{k})$ is a modified resolvent operator. Its singular
value decomposition can be written as

$\mathcal{H}^{e}(\bm{k})=\hat{\bm{\Psi}}^{e}(\bm{k})\bm{\Sigma}^{e}(\bm{k})\hat{\bm{\Phi}}^{*,e}(\bm{k}).$
(12)

The superscript $e$ differentiates the eddy resolvent modes, henceforth
referred to as eddy modes, and singular values from their standard resolvent
counterparts in Eq. (6). The interpretation of each term in the decomposition
is similar to Eq. (6) in that $\hat{\Psi}^{e}(\bm{k})$ consists of orthogonal
basis functions for the velocity field and the diagonal matrix
$\bm{\Sigma}^{e}(\bm{k})$ ranks the $p$th structure by its gain using an inner
product proportional to its kinetic energy. The matrix
$\hat{\Phi}^{e}(\bm{k})$, on the other hand, contains orthogonal basis
functions for the disturbance field $\hat{\bm{d}}(\bm{k})$ which is less
interpretable than $\hat{\bm{f}}(\bm{k})$. Despite this drawback, the addition
of eddy viscosity is expected to partially model the effect of
$\hat{\bm{f}}(\bm{k})$ and thus improve the efficiency of eddy modes as a
basis for the velocity field.

###  II.4 Spectral proper orthogonal decomposition

The efficiency of resolvent and eddy modes as a basis for the velocity field
can be assessed by projecting them onto SPOD modes, which are computed
directly from data. The SPOD modes are computed with the same procedure
described in Ref. Towne  _et al._ (2018) so only a brief summary is presented
here. Using Welch’s method (Welch, 1967), the DNS data for a particular
$(k_{x},k_{y})$ are divided into overlapping segments containing 512 snapshots
with 75% overlap. Each segment is Fourier-transformed in time and the Fourier
modes for a specific frequency $\omega$ can be arranged into the new data
matrix

$\hat{\bm{Q}}(\bm{k})=\left[\begin{array}[]{cccc}\hat{\bm{q}}_{\omega}^{(1)}&\hat{\bm{q}}_{\omega}^{(2)}&\cdots&\hat{\bm{q}}_{\omega}^{(s)}\end{array}\right]\in\mathbb{C}^{m\times
s},$ (13)

where $m$ represents the number of states and $s$ the number of segements. The
cross-spectral density matrix for a specific wavenumber triplet
$\hat{\bm{S}}(\bm{k})$ is

$\hat{\bm{S}}(\bm{k})=\hat{\bm{Q}}(\bm{k})\hat{\bm{Q}}^{*}(\bm{k}).$ (14)

The SPOD eigenvectors (or modes) $\hat{\bm{V}}(\bm{k})$ and eigenvalues
$\bm{\Lambda}(\bm{k})$ can be obtained by performing an eigenvalue
decomposition of the cross-spectral density matrix

$\hat{\bm{S}}(\bm{k})\hat{\bm{V}}(\bm{k})=\hat{\bm{V}}(\bm{k})\bm{\Lambda}(\bm{k}).$
(15)

##  III DNS dataset

$Re_{\tau}$ | $L_{x}$ | $L_{y}$ | $L_{z}$ | $N_{x}$ | $N_{y}$ | $N_{z}$ | $\Delta x^{+}$ | $\Delta y^{+}$ | $\Delta z^{+}_{min}$ | $\Delta z^{+}_{max}$
---|---|---|---|---|---|---|---|---|---|---
550 | $2\pi$ | $\pi$ | $2h$ | 256 | 256 | 201 | 13.5 | 6.75 | $6.79\times 10^{-2}$ | 8.64
Table 1: Channel flow DNS parameters.

A DNS of channel flow at $Re_{\tau}=550$ is performed using the ChannelFlow
pseudo-spectral code Gibson  _et al._ (2019). Table 1 summarizes the
parameters of the simulation which was solved on a domain with dimensions
$2\pi\times\pi\times 2h$ in the streamwise ($L_{x}$), spanwise ($L_{y}$) and
wall-normal ($L_{z}$) directions. There are $N_{x}=N_{y}=256$ equally spaced
points in the streamwise and spanwise directions and $N_{y}=201$ points in the
wall-normal direction on a Chebyshev grid. Periodic boundary conditions are
employed in the streamwise and spanwise directions while no-slip boundary
conditions are enforced on the channel walls. Further details on the mesh
discretization $(\Delta x^{+},\Delta y^{+},\Delta z^{+}_{min},\Delta
z^{+}_{\max})$ are presented in Table 1. The mean velocity and Reynolds stress
profiles are presented in Figs. 1(a) and (b), respectively. All profiles show
good agreement with the DNS results from Ref. Lee and Moser (2015) despite the
smaller computational box in this study.

Figure 1: (a) Mean velocity profile and (b) Reynolds stress profiles from the
present study (solid lines) and Ref. Lee and Moser (2015) (open squares).

SPOD is performed on a database of 6784 snapshots at a time resolution of
$\Delta t=0.2$. The data are divided into equal segments containing 512
snapshots with an overlap of 75% resulting in 50 blocks. The cross-spectral
density matrices in Eq. (14) are computed using Welch’s method with a Hamming
window. The SPOD modes and their respective energies for a desired frequency
are obtained from the eigenvectors and eigenvalues, respectively, of the
cross-spectral density matrices.

##  IV Resolvent and eddy analysis predictions

In this section, the predictions from resolvent and eddy analysis, i.e. the
resolvent supplemented with an eddy viscosity model, are compared to DNS data
using several measures. To begin with, the low-rank behavior of the resolvent
is compared to the turbulent kinetic energy spectrum of DNS in Sec. IV.1.
Next, the projection of the leading resolvent and eddy modes onto the dominant
SPOD mode is computed in Sec. IV.2 for the most energetic wavenumber pairs.
Structures convecting at a range of wave speeds $c^{+}=\omega/k_{x}$ are
considered, and projections are computed independently for each wave speed. In
Sec. IV.3, the most energetic wave speed for a range of $(k_{x},k_{y})$
wavenumber pairs in DNS is compared to the most amplified wave speed
identified by resolvent and eddy analysis. This motivates the selection of
specific wavenumber triplets to examine in greater detail in Sec. V in which
the SPOD, resolvent, and eddy mode shapes are compared directly.

###  IV.1 Low-rank maps

An important aspect of the resolvent operator is its rank or, more
specifically, the ratio of the total energy that is captured by the leading
resolvent modes alone. Reference Moarref  _et al._ (2013) was able to show
that for a turbulent channel flow at $Re_{\tau}=2003$, in the region of the
$(k_{x},k_{y})$ space where the real flow is energetic (as seen from
experiments or DNS), the standard resolvent operator tends to be low-rank.
With this observation in mind, in this section the low-rank maps from
resolvent and eddy analysis are compared to the turbulent kinetic energy
spectra from DNS (data provided by Refs. del Álamo and Jiménez (2003) and del
Álamo  _et al._ (2004)). As noted by Ref. Moarref  _et al._ (2013), the
symmetry of the channel leads to the resolvent singular values coming in
(approximately) equal pairs. One mode in each pair is symmetric with respect
to the channel centerline while the other mode is anti-symmetric. The low-rank
behavior of the resolvent for the channel can therefore be studied by
computing the ratio of the dominant pair of singular values to the sum of all
singular values, i.e.

$\mathcal{R}(\bm{k})=\frac{\sigma_{1}^{2}(\bm{k})+\sigma_{2}^{2}(\bm{k})}{\sum_{p}\sigma_{p}^{2}(\bm{k})}.$
(16)

Since $\mathcal{R}$ depends on a wavenumber triplet, $\mathcal{R}$ is computed
across $(\lambda_{x},\lambda_{y})$ for a fixed wave speed in order to
facilitate visualization.

Figure 2: Low-rank maps for (a) resolvent and (b) eddy analysis for a fixed
wave speed of $c^{+}=10$. Contours of the turbulent kinetic energy spectrum at
$z^{+}=15$ from Refs. del Álamo and Jiménez (2003) and del Álamo  _et al._
(2004) are denoted in black.

The low-rank maps for $c^{+}=10$ are presented in Figs. 2(a) and (b) for
resolvent and eddy analysis, respectively. The wall-normal location
corresponding to $U^{+}=10$ is approximately $z^{+}=15$ as indicated by the
dash-dotted line in Fig. 1(a). The colours in the figure represent the value
of $\mathcal{R}(\bm{k})$, while the black contour lines represent the
streamwise turbulent kinetic energy spectrum from DNS at $z^{+}\approx 15$ del
Álamo and Jiménez (2003); del Álamo  _et al._ (2004). Figure 2(a) shows that
there is very good agreement between the scales with the most energy in DNS
and the scales at which the resolvent is low-rank. This is to be expected
since the energetic structures in the flow tend to arise from linear
amplification mechanisms that are identified well by resolvent analysis
Jovanović and Bamieh (2005); McKeon and Sharma (2010). The eddy low-rank map
in Fig. 2(b) has two main differences from its resolvent counterpart in Fig.
2(a). The first is that there is less agreement between the low-rank map and
the energy spectra. The most energetic scales from DNS with spanwise
wavelengths of $\lambda_{y}^{+}\approx 100$, for example, do not coincide with
the peak in the eddy low-rank map which occurs at $\lambda_{y}^{+}\approx 80$.
The second is that there are two peaks in the low-rank map for the eddy
analysis whereas there is only one for resolvent analysis. The two peaks from
the eddy operator occur at $\lambda_{y}^{+}=80$ and $\lambda_{y}=3.5h$,
wavelengths that have been previously observed by Refs. del Álamo and Jiménez
(2006) and Hwang and Cossu (2010) in the context of transient growth and
harmonic forcing analyses, respectively. It can be observed that there is some
energy from larger scales present at this wall-normal location as attested to
by the closed contour of kinetic energy centered on
$(\lambda_{x},\lambda_{y})=(14h,2h)$ or
$(\lambda_{x}^{+},\lambda_{y}^{+})=(7700,1100)$. Nevertheless, the Reynolds
number of this study is too low to conclude that the second peak around
$\lambda_{y}=3.5h$ is indicative of the footprint of large-scales at the wall.

Figure 3: Low-rank maps for (a) resolvent and (b) eddy analysis for a fixed
wave speed of $c^{+}=18.9$. Contours of the turbulent kinetic energy spectrum
at $z^{+}=-h/2$ from Refs. del Álamo and Jiménez (2003) and del Álamo  _et
al._ (2004) are denoted in black.

The low-rank maps for $c^{+}=18.9$ are shown in Fig. 3. Similar to the
previous wave speed, there is good agreement between the energy spectrum from
DNS and the low-rank map for resolvent analysis in Fig. 3(a). The agreement
between the DNS and low-rank map for eddy analysis in Fig. 3(b), on the other
hand, is worse than it was for the previous wave speed. The only area of
slight agreement occurs for scales with the longest streamwise wavelengths
from DNS. Perhaps the most remarkable aspect of the low-rank map in Fig. 3(b)
is that it closely resembles the low-rank map in Fig. 2(b). The resolvent low-
rank maps, on the other hand, are influenced significantly by the wave speed.
It can be concluded that the wave speed does not have a significant impact on
the eddy low-rank maps. For both wave speeds considered, there are two
spanwise wavelengths for which the eddy operator is low-rank. The only
difference is that $\mathcal{R}$ for $c^{+}=18.9$ decreases for the smaller
spanwise wavelength and increases for the larger spanwise wavelength relative
to $\mathcal{R}$ for $c^{+}=10$.

It can be remarked that the addition of eddy viscosity significantly distorts
the linear mechanisms identified by resolvent analysis. As such, eddy analysis
is less successful in identifying the energetic scales for a specified wave
speed. These results, notwithstanding, do not quantify the accuracy of
resolvent or eddy analysis in predicting flow structures in turbulent channel
flow. It will be seen that the two peak spanwise wavelengths identified by the
low-rank maps provide a valuable clue in identifying the types of structures
that eddy analysis predicts with good accuracy.

###  IV.2 Projection of resolvent modes onto SPOD modes

The objective of this section is to quantify the accuracy of resolvent and
eddy analysis by projecting the leading SPOD mode $\hat{\bm{v}}_{1}(\bm{k})$
from DNS onto the leading resolvent $\hat{\bm{\psi}}_{1}(\bm{k})$ and eddy
$\hat{\bm{\psi}}^{e}_{1}(\bm{k})$ modes. Similar analyses have been performed
by Refs. Abreu  _et al._ (2020) and Pickering  _et al._ (2021) for turbulent
pipe flow and turbulent jets, respectively. To account for the pairing of
resolvent and eddy modes, the leading SPOD mode is projected onto both the
first and second resolvent/eddy modes

$\gamma(\bm{k})=\sqrt{\left(\frac{\left<\hat{\bm{v}}_{1}(\bm{k}),\hat{\bm{\psi}}_{1}(\bm{k})\right>}{\|\hat{\bm{v}}_{1}(\bm{k})\|\cdot\|\hat{\bm{\psi}}_{1}(\bm{k})\|}\right)^{2}+\left(\frac{\left<\hat{\bm{v}}_{1}(\bm{k}),\hat{\bm{\psi}}_{2}(\bm{k})\right>}{\|\hat{\bm{v}}_{1}(\bm{k})\|\cdot\|\hat{\bm{\psi}}_{2}(\bm{k})\|}\right)^{2}},$
(17a)
$\gamma^{e}(\bm{k})=\sqrt{\left(\frac{\left<\hat{\bm{v}}_{1}(\bm{k}),\hat{\bm{\psi}}^{e}_{1}(\bm{k})\right>}{\|\hat{\bm{v}}_{1}(\bm{k})\|\cdot\|\hat{\bm{\psi}}^{e}_{1}(\bm{k})\|}\right)^{2}+\left(\frac{\left<\hat{\bm{v}}_{1}(\bm{k}),\hat{\bm{\psi}}^{e}_{2}(\bm{k})\right>}{\|\hat{\bm{v}}_{1}(\bm{k})\|\cdot\|\hat{\bm{\psi}}^{e}_{2}(\bm{k})\|}\right)^{2}}.$
(17b)

Both projection coefficients $\gamma(\bm{k})$ and $\gamma^{e}(\bm{k})$ have a
maximum value of unity, which indicates perfect alignment between SPOD and
resolvent/eddy modes. A value of zero indicates that the mode shapes are
orthogonal.

Figure 4: Projection of the leading resolvent mode onto the leading SPOD mode
for the wave speeds $c^{+}=10.3,~{}12.0,~{}13.7,~{}15.5,~{}17.2,~{}18.9$ in
ascending order. Figure 5: Projection of the leading eddy mode onto the
leading SPOD mode for the wave speeds
$c^{+}=10.3,~{}12.0,~{}13.7,~{}15.5,~{}17.2,~{}18.9$ in ascending order.

Figure 4 illustrates $\gamma$ for wavenumbers that satisfy $1\leq k_{x}\leq
19$ and $0\leq k_{y}\leq 40$. These wavenumber pairs contain most of the
kinetic energy in the flow and include structures associated with the near-
wall cycle. Streamwise-constant modes are excluded since the wave speed is
undefined for $k_{x}=0$. They will be discussed in greater detail in Sec. V.
The colorbar in Fig. 4 is restricted to a range of $[0.5,~{}1]$ to facilitate
identification of wavenumber pairs where there is significant overlap between
the SPOD and resolvent modes. Each panel in Fig. 4 represents a different wave
speed. The lowest wave speed considered is $c^{+}=10.3$ in Fig. 4(a) and
$c^{+}$ increases at increments of approximately 1.7, culminating in a wave
speed of $c^{+}=18.9$ in Fig. 4(f). The lower bound on $c^{+}$ is motivated by
the near-wall streaks being most energetic at a wall-normal location of
$z^{+}=15$ where $U^{+}\approx 10$ McKeon and Sharma (2010). The upper bound
on $c^{+}$ is chosen since it is approximately $c^{+}=U^{+}_{CL}-2$, which
translates to a wall-normal location of $z=-h/2$. Past studies Morra  _et al._
(2019); Symon  _et al._ (2020); Morra  _et al._ (2021) have noted good
agreement between SPOD and eddy modes for structures that are most energetic
at this wall-normal location. The eddy viscosity profile $\nu_{T}(z)$ also
reaches a maximum at $z=-h/2$, an observation that is shown to be significant
in Sec. V.

Figure 4 shows that there is good agreement between the leading SPOD and
resolvent modes for $c^{+}=10.3$. Despite some outliers, Fig. 4(a) indicates
that the highest projections are for modes that satisfy $k_{x}<k_{y}$. This is
consistent with the results of Ref. Abreu  _et al._ (2020) who noted that the
lift-up mechanism leads to large amplification for high aspect ratio scales
where $\AR=k_{y}/k_{x}$. As the wave speed increases, however, the values of
$\gamma$ decline quite significantly, indicating poorer predictions of the
mode shapes from resolvent analysis. Figure 5 presents the projection
coefficient for eddy analysis. It is observed that there is significantly
better agreement between the leading SPOD and eddy modes at all wave speeds in
comparison to the standard resolvent in Fig. 4. There are, however, some
similarities between $\gamma$ and $\gamma^{e}$. First, the highest values of
$\gamma^{e}$ are obtained for high aspect ratio structures. Second, the wave
speed has a major influence on $\gamma^{e}$. In general, the SPOD and eddy
modes overlap less as the wave speed is increased. Unlike $\gamma$, however,
there are several wavenumber pairs for which $\gamma^{e}>0.9$ at the largest
wave speed. The largest $\gamma^{e}$ are clustered around the largest scales
in Fig. 5(f) whereas in Fig. 5(a), they are centered around the near-wall
cycle mode $(k_{x},k_{y})=(4,30)$.

###  IV.3 Frequency response

The previous sections show that there is a trade-off when eddy viscosity is
added to the resolvent operator. On one hand, eddy viscosity distorts the
linear dynamics of the operator such that the low-rank map resembles less the
turbulent kinetic energy spectrum. On the other hand, the eddy modes have
larger projections onto SPOD modes than their resolvent mode counterparts.
Although these comparisons have been studied for a variety of wave speeds, the
impact of eddy viscosity on the resolvent frequency response has yet to be
analyzed. In order to be consistent with previous sections, the frequency
response is considered from a wave speed point of view. The impact of eddy
viscosity can be better appreciated by considering a single wavenumber pair
initially before investigating all energetic wavenumber pairs as done in the
previous sections. The wavenumber pair $(k_{x},k_{y})=(4,30)$, which
corresponds to $(\lambda_{x}^{+},\lambda_{y}^{+})=(864,116)$, is selected as
it corresponds to the near-wall streaks in the DNS used in this study.

To gain a better understanding of what the addition of eddy viscosity does to
modify the frequency response of the resolvent operator, we look at a range of
eddy viscosity based models, where the strength of the eddy viscosity is
gradually increased from zero (equivalent to the standard resolvent) to the
full eddy-viscosity. The strength of the eddy viscosity is adjusted
artificially by introducing the scaling factor $\mathcal{S}$ such that Eq.
(10) becomes

$\nu_{T}(z)=\frac{\nu}{2}\left(1+\mathcal{S}\left[\frac{\kappa}{3}(1-z^{2})(1+2z^{2})(1-\text{exp}(|z-1|Re_{\tau}/A)\right]^{2}\right)^{1/2}+\frac{\nu}{2},$
(18)

where $\mathcal{S}\in[0,~{}1]$. Setting $\mathcal{S}=0$ or $\mathcal{S}=1$ is
equivalent to resolvent analysis or eddy analysis, respectively. A similar
parameter was introduced by Ref. Gupta  _et al._ (2021) to derive a scale-
dependent eddy viscosity for linear estimation of a turbulent channel flow at
$Re_{\tau}=2003$.

Figure 6: Most energetic wave speed for $(k_{x},k_{y})=(4,30)$ as a function
of different strengths of the eddy viscosity as denoted by $\mathcal{S}$.

The first singular value $\sigma_{1}$ is plotted against $c^{+}$ in Fig. 6 for
various strengths of eddy viscosity, i.e. different values of $\mathcal{S}$.
As the value of $\mathcal{S}$ increases, both the maximum amplification and
the most amplified wave speed decrease. For this particular scale, moreover,
the maximum amplification declines by a factor of 10 and the most amplified
wave speed slows down substantially from $c^{+}=15.5$ to $c^{+}=9.5$. A crude
explanation for this behavior is that the damping supplied by eddy viscosity
results in slower, less amplified structures. Since resolvent analysis has
absolutely no damping other than molecular viscosity, the structures are
allowed to convect more quickly. It can, therefore, be expected that for an
arbitrary scale the most amplified wave speed predicted by resolvent analysis
is going to be greater than that predicted by eddy analysis.

Figure 7: Most energetic wave speed computed from (a) DNS compared to the most
amplified wave speed predicted by (b) resolvent analysis and (c) eddy
analysis.

This hypothesis is tested in Fig. 7, which compares the most energetic wave
speed from DNS to the most amplified wave speed predicted by resolvent and
eddy analysis. The only similarity among the three panels is that the most
energetic/amplified wave speed is primarily governed by the spanwise
wavenumber. Wider structures, i.e. those with small spanwise wavenumbers,
travel faster than relatively less wide structures. The streamwise wavenumber
plays a bigger role in DNS and resolvent analysis than it does for eddy
analysis. Figures 7(a) and (b) show that for fixed $k_{y}$, the most
energetic/amplified wave speed increases as a function of $k_{x}$. Thus, the
trend for $k_{x}$ is different from $k_{y}$ in that longer structures, i.e.
those with smaller $k_{x}$, travel slower, while wider structures, i.e. those
with smaller $k_{y}$ travel faster. Another interpretation of these trends is
that higher aspect ratio structures convect more slowly and are thus more
energetic closer to the wall. Lower aspect ratio structures, meanwhile,
convect more quickly and are thus more energetic away from the wall. These
trends are consistent with observations from Ref. Hwang (2015) among others
who report aspect ratios of approximately $\AR=8$ for near-wall coherent
motions and $\AR=2-3$ for large-scale outer motions.

Figure 8: Difference between the most amplified wave speed predicted by (a)
resolvent and (b) eddy analysis compared to the most energetic wave speed
computed in DNS.

Figure 8 plots the difference between the most amplified wave speed predicted
by resolvent/eddy analysis and the most energetic wave speed computed from
DNS. Red squares indicate that the most amplified wave speed is too high and
blue squares that the most amplified wave speed is too low. It is striking how
large the differences are particularly for the relatively small Reynolds
number considered in this study which restricts the range of wave speeds that
can be expected in the flow. Figure 8(a) shows that the most amplified
structures in resolvent analysis travel faster than their true speeds whereas
in Fig. 8(b), the most amplified structures in eddy analysis travel slower.
Despite the large differences between DNS and the linear analyses, there are
two regions in Fig. 8(b) where the discrepancy is small. The first is denoted
by the dashed, black rectangle and corresponds to small values of both $k_{x}$
and $k_{y}$, where $c^{+}\approx 18.9$. The second is denoted by the solid
black rectangle, which contains wavenumber pairs in the range
$(k_{x},k_{y})=(1\leq k_{x}\leq 5,30\leq k_{y}\leq 34)$ where $c^{+}\approx
10$. As will be discussed in the next section, the reasons for why these two
wavenumber regions are estimated reasonably well by the eddy model can be
related back to the wall-normal profile of the eddy viscosity $\nu_{T}(z)$.

##  V Mode shapes and the Cess eddy viscosity profile

The previous section analyzed the linear predictions of resolvent and eddy
analysis using scalar quantities such as the low-rank maps, projection
coefficients, and most amplified wave speeds. This section considers the wall-
normal profiles of both mode shapes in Sec. V.1 and the Cess eddy viscosity
profile itself in Sec. V.2 to predict the agreement between SPOD and
resolvent/eddy modes at higher Reynolds numbers.

###  V.1 SPOD, resolvent, and eddy modes

Mode | $k_{x}$ | $k_{y}$ | $c^{+}_{\text{max,SPOD}}$ | $c^{+}_{\text{max,resolvent}}$ | $c^{+}_{\text{max,eddy}}$ | $\gamma$ | $\gamma^{e}$ | $\gamma_{95}$ | $\gamma_{95}^{e}$
---|---|---|---|---|---|---|---|---|---
0 | 0 | 4 | $\omega=0$ | $\omega=0$ | $\omega=0$ | 0.898 | 0.978 | 6 | 2
1 | 1 | 2 | 18.5 | 19.0 | 18.5 | 0.663 | 0.956 | 36 | 2
2 | 2 | 8 | 16.8 | 20.5 | 13.8 | 0.586 | 0.907 | 56 | 8
3 | 4 | 30 | 10.1 | 15.3 | 9.50 | 0.864 | 0.959 | 40 | 2
Table 2: The wavenumber triplets for modes 0-3 along with predictions from
resolvent and eddy analysis.

The SPOD, resolvent, and eddy mode shapes are compared for four wavenumber
triplets that are described in Table 2. Mode 0 is a streamwise-constant mode
and its most energetic frequency is $\omega=0$. The other wavenumber pairs are
chosen such that Modes 1, 2, and 3 have their peak energy around $z/h=-h/2$,
$z/h=-h/4$, and $z^{+}=15$, respectively. The most energetic wave speed from
SPOD is reported as $c_{\text{max,SPOD}}^{+}$. Table 2 reports
$c^{+}_{\text{max,resolvent}}$ and $c^{+}_{\text{max,eddy}}$, the most
amplified wave speeds identified by resolvent and eddy analysis, respectively.
The wave speed selected for plotting purposes, however, is kept fixed at
$c_{\text{max,SPOD}}^{+}$. The purpose of comparing $c_{\text{max,SPOD}}^{+}$
and $c^{+}_{\text{max,eddy}}$ is to emphasize that eddy analysis produces good
predictions when $c_{\text{max,SPOD}}^{+}\approx c^{+}_{\text{max,eddy}}$. The
parameter $c_{\text{max,resolvent}}^{+}$, meanwhile, does not provide
information about the quality of predictions from resolvent analysis. Finally,
Table 2 quantifies the number of resolvent/eddy modes needed such that the
velocity from DNS is reconstructed to $95\%$ accuracy. In other words, the
number of modes required to obtain a value of $\gamma=0.95$ using Eq. 17.

Figure 9: The (a) streamwise, (b) spanwise, and (c) wall-normal component of
the leading SPOD, resolvent, and eddy mode for Mode 0 for which
$\bm{k}=(0,4,\omega=0)$.

Figure 9 presents the SPOD, resolvent, and eddy mode for Mode 0. As indicated
in Table 2, both resolvent and eddy analysis correctly identify $\omega=0$ as
the most amplified frequency. The projection coefficients $\gamma$ and
$\gamma^{e}$ are also large for this choice of $\bm{k}$. The number of
resolvent or eddy modes needed to achieve $\gamma=0.95$, therefore, is
$\gamma_{95}=6$ and $\gamma^{e}_{95}=2$, respectively. The agreement between
the resolvent and SPOD modes themselves, however, is not as compelling as
$\gamma=0.898$ might suggest. Although the streamwise velocity component is
predicted reasonably well, the spanwise and wall-normal components are
significantly underestimated due to the high non-normality of the resolvent
operator Trefethen  _et al._ (1993); Schmid and Henningson (2001). The eddy
operator is also non-normal but the addition of eddy viscosity results in a
more normal operator (Symon  _et al._ , 2021). Since the streamwise component
is also dominant for the SPOD mode, it has a disproportionate influence on
$\gamma$ thus resulting in a higher value than might be expected from visual
inspection of the mode shapes. The agreement between the SPOD and eddy modes,
on the other hand, is good for all velocity components, resulting in
$\gamma^{e}=0.978$.

Figure 10: The (a) streamwise, (b) spanwise, and (c) wall-normal component of
the leading SPOD, resolvent, and eddy mode for Mode 1 for which
$\bm{k}=(1,2,18.5)$.

Mode 1 is an energetic structure in the outer region of the flow with a
wavenumber triplet of $\bm{k}=(1,2,18.5)$. Table 2 shows that
$c^{+}_{\text{max,SPOD}}=18.5$ and this is in good agreement with predictions
from resolvent and eddy analysis. Figure 10 illustrates that while the
resolvent modes do not agree well with the SPOD modes, the eddy modes show
reasonable agreement with SPOD. This is reflected in the projection
coefficients as $\gamma^{e}=0.956$ is higher than $\gamma=0.663$. It also
takes significantly fewer eddy modes to reconstruct the leading SPOD mode
since as $\gamma^{e}_{95}=2$ whereas $\gamma_{95}=36$. Refs. Rosenberg and
McKeon (2019); Morra  _et al._ (2019); Symon  _et al._ (2021) also observed
that, when using resolvent modes, many suboptimal modes are needed to
reconstruct the velocity field for high aspect ratio structures. One reason is
that the streamwise and spanwise components of velocity are highly localized
in the wall-normal direction due to the critical-layer mechanism McKeon and
Sharma (2010); Morra  _et al._ (2019); Madhusudanan  _et al._ (2019); Vadarevu
_et al._ (2019); Symon  _et al._ (2020). The eddy modes, on the other hand,
are smoothed out in $z$ by the eddy viscosity. Another factor is that, similar
to Mode 0, the wall-normal velocity component is underestimated by resolvent
analysis as seen in Fig. 10(c). The agreement between the SPOD and eddy modes,
meanwhile, is very good for all velocity components and wall-normal locations
except the near-wall region of the streamwise component, as seen in Fig.
10(a).

Figure 11: The (a) streamwise, (b) spanwise, and (c) wall-normal component of
the leading SPOD, resolvent, and eddy mode for Mode 2 for which
$\bm{k}=(2,8,16.8)$.

Mode 2 is another structure in the outer region of the flow with a wavenumber
triplet of $\bm{k}=(2,8,16.8)$. The streamwise velocity component is most
energetic at $z/h=-0.75$ which is closer to the wall in comparison to
$z/h=-0.5$ for Mode 1. Neither resolvent nor eddy analysis are capable of
predicting the most energetic wave speed for Mode 2. Consistent with the
trends observed in Fig. 8, resolvent analysis predicts a wave speed that is
too fast while eddy analysis predicts a wave speed that is too slow. As
explained earlier, for plotting the resolvent and eddy modes in Fig. 11,
$c^{+}_{resolvent}$ and $c^{+}_{eddy}$ are chosen to be the maximum wave-speed
as identified by SPOD, i.e. $c^{+}_{max,SPOD}$. It can be remarked that the
projection coefficients belie the true agreement between various mode shapes.
When considering the resolvent modes, although $\gamma=0.586$ seems low, the
modes still capture the wall-normal location of the peak streamwise energy
reasonably well. However, similar to Mode 1, the modes are localized about the
critical-layer and the wall-normal component is underestimated resulting in a
relatively low $\gamma$. The eddy modes appear to be in significantly better
agreement as $\gamma^{e}=0.907$. Unlike resolvent analysis, however, the wall-
normal location of the peak in streamwise and spanwise energies of the eddy
modes fall below their true locations, a trend that becomes more apparent for
higher Reynolds numbers (see Ref. Symon  _et al._ (2020)). Even though it is
not explicitly shown here for the sake of brevity, for structures with
$c^{+}_{\text{max,eddy}}<c^{+}_{\text{max,SPOD}}$, the peak streamwise and
spanwise energies of the eddy mode are located below the correct wall-normal
location identified by SPOD. Mode 2, consequently, is the only wavenumber
triplet for which $\gamma^{e}_{95}=8$ since suboptimal modes are required to
“lift” the structure to its proper height. 58 modes, meanwhile, are required
for resolvent modes to achieve $\gamma=0.95$.

Figure 12: The (a) streamwise, (b) spanwise, and (c) wall-normal component of
the leading SPOD, resolvent, and eddy mode for Mode 3 for which
$\bm{k}=(4,30,10.1)$.

Figure 12 considers Mode 3, which is representative of the near-wall cycle.
The agreement between the SPOD and resolvent wave speeds is poor but agreement
between the mode shapes is high as $\gamma=0.864$. Similar to Mode 0, the
streamwise velocity component is the dominant component for Mode 3, thus
relegating the influence of the other two velocity components in computing
$\gamma$. As seen in Table 2, eddy analysis is close to identifying the most
energetic wave speed. Furthermore, the SPOD and eddy mode shapes in Fig. 12,
are in excellent agreement for the streamwise velocity component and good
agreement for the spanwise and wall-normal velocity components, resulting in
$\gamma^{e}=0.959$. It is clear, furthermore, that eddy modes are a more
efficient basis since $\gamma^{e}_{95}=2$ compared to $\gamma_{95}=40$.
Similar agreement between DNS and eddy analysis predictions for near-wall
structures has also been reported by Refs. Morra  _et al._ (2019, 2021).

###  V.2 Cess eddy viscosity and effective Reynolds number

Figure 13: Wall-normal profiles of (a) $Re_{\tau}\nu_{T}$ and (b) $\nu_{T}$
for various $Re_{\tau}$.

In this section, the Cess eddy viscosity profile is studied in greater detail
in order to predict the agreement between SPOD and resolvent/eddy modes at
higher Reynolds numbers. The profiles are obtained for a range of Reynolds
numbers using Eq. (10). It is particularly interesting to look at the profiles
around $z^{+}=15$ and $z=-h/2$ due to the good agreement between SPOD and eddy
analysis for Modes 1 and 3 that have their peak energies at these wall-
heights. In Fig. 13(a), the eddy viscosity profiles are pre-multiplied by
$Re_{\tau}$, resulting in a collapse near the wall. The profiles for
$Re_{\tau}=180$ and $Re_{\tau}=550$ begin to diverge from the other Reynolds
numbers around $z^{+}=10$ and $z^{+}=20$, respectively, but the agreement
among all profiles is good, particularly around $z^{+}=15$. This wall-normal
location is important since it coincides with the location of the structures
associated with the near-wall cycle that travel at a wave speed of $c^{+}=10$.
It therefore seems reasonable to assume that, for all Reynolds numbers, the
eddy analysis will correctly predict the structures that are located at these
wall-heights of $z^{+}\approx 15$ and therefore convect at wave speeds around
$c^{+}=10$. Consistent with this observation, Ref. Morra  _et al._ (2019)
observed that at $Re_{\tau}=1007$, the eddy analysis is able to get reasonable
predictions for the structure $(k_{x},k_{y})=(14,63)$ convecting at
$c^{+}=10$.

In Fig. 13(b), $\nu_{T}(z)$ is plotted in outer units for Reynolds numbers in
the range $180\leq Re_{\tau}\leq 20000$. The profiles for all Reynolds numbers
other than $Re_{\tau}=180$ are virtually indistinguishable in the outer
region. The largest differences occur near and at the walls where
$\nu_{T}(0)=\nu_{T}(2h)=\nu$. The maximum value of $\nu_{T}$ is 0.08 and
occurs at $z=\pm h/2$. Mode 1, for which the eddy analysis gives good
predictions, was also found to be most energetic at this wall-height of $z=\pm
h/2$. Since $\nu_{T}$ at $z=\pm h/2$ remains roughly constant with Reynolds
number, it can be hypothesized that the eddy mode shapes for wavenumbers with
$c^{+}_{\text{max}}=U^{+}(\pm h/2)\approx U^{+}_{CL}-2$ are unaffected by
$Re_{\tau}$. This hypothesis is tested in Fig. 14 where the mode shapes of
structures with $c^{+}_{\text{max}}\approx U^{+}_{CL}-2$ are compared across
different $Re_{\tau}$. A second hypothesis, which is also tested in Fig. 14,
is that the wall-normal-varying eddy viscosity $\nu_{T}(z)$ can be replaced by
a constant eddy viscosity model, i.e. $\nu_{T}=Re_{T}$, where $Re_{T}$ is the
effective Reynolds number. For the case of a turbulent jet, References
Pickering  _et al._ (2021); Kuhn  _et al._ (2021, 2022) showed that a constant
eddy viscosity model can improve the agreement between SPOD and eddy modes.

Figure 14: Mode shapes for $\bm{k}=(1,2,U^{+}_{CL}-2)$ from (a) resolvent
analysis, (b) eddy analysis, and (c) eddy analysis setting $Re_{T}=12.5$. The
Reynolds numbers range from $Re_{\tau}=180$ up to $Re_{\tau}=20000$.

Figure 14 compares the mode shapes for $\bm{k}=(1,2,U^{+}_{CL}-2)$ from
resolvent analysis, eddy analysis, and eddy analysis setting
$Re_{T}=1/(\text{max}(\nu_{T}))=12.5$ as done in Hwang (2016). If
$\nu_{T}=Re_{T}$, then only the mean profile is affected by changes in
Reynolds number. All sets of modes are computed for the same Reynolds numbers
that appeared in Fig. 13 and are normalized by the maximum value of the
streamwise velocity component. The resolvent modes in Fig. 14(a) are
influenced heavily by the choice of $Re_{\tau}$. The streamwise and spanwise
velocity components become increasingly localized about the critical layer at
$z=-h/2$ as the Reynolds number increases. The shape of the wall-normal
component, on the other hand, is roughly constant but its magnitude relative
to the wall-parallel velocity components decreases with increasing Reynolds
number.

The impact of Reynolds number on the eddy and constant $Re_{T}$ modes in Figs.
14(b,c) is negligible. The only difference among the eddy modes in Fig. 14(b)
is that for $Re_{\tau}=180$, the streamwise and spanwise velocity components
are less attached to the wall in comparison to the other profiles which appear
more blunt, i.e. flatter, near the wall. A similar difference emerges between
the eddy and constant $Re_{T}$ modes in that neither the streamwise nor
spanwise components exhibit blunt behavior near the wall for any Reynolds
number considered. The agreement between eddy and constant $Re_{T}$ modes,
nonetheless, is remarkable given that the constant eddy viscosity model is so
simple.

The applicability of the constant eddy viscosity model, however, is limited to
scales that have maximum streamwise energy at $z/h=\pm 0.5$. Figure 15
compares the mode shapes for
$(k_{x}^{+},k_{y}^{+},c^{+})=(2\pi/1000,2\pi/100,10)$, which are
representative of the near-wall cycle, from resolvent analysis, eddy analysis,
and eddy analysis setting $Re_{T}=1/(\text{max}(\nu_{T}))=12.5$. The Reynolds
number has little impact on the resolvent and eddy modes but has a major
impact on the constant $Re_{T}$ modes. As the Reynolds number increases, the
mode shapes become increasingly less localized about the critical layer. The
location of the peak energy of all three velocity components, furthermore,
gradually shifts closer to the channel centerline and away from the wall.
Although it is not shown in the interest of brevity, the same trends can be
observed for other values of $Re_{T}$. The wall-varying profile of the eddy
viscosity near the wall, consequently, is essential to capture the correct
mode shapes for wavenumber triplets that are associated with the near-wall
cycle.

Figure 15: Mode shapes for
$(k_{x}^{+},k_{y}^{+},c^{+})=(2\pi/1000,2\pi/100,10)$ from (a) resolvent
analysis, (b) eddy analysis, and (c) eddy analysis setting $Re_{T}=12.5$. The
Reynolds numbers range from $Re_{\tau}=180$ up to $Re_{\tau}=20000$.

##  VI Nonlinear Energy transfer and eddy viscosity

The improved agreement between SPOD and eddy modes indicates that eddy
viscosity is able to partially model the effect of $\hat{\bm{f}}$ for
particular scales such as Modes 1 and 3. In this section, the energy transfers
that are introduced by the eddy viscosity, herein referred to as eddy
dissipation, are examined in greater detail to explain the success or failure
of eddy analysis in predicting the correct structures. Section VI.1 introduces
eddy dissipation which consists of two terms. The first originates from the
wall-normal-varying effective Reynolds number $\nu_{T}(z)-\nu$ and the second
from the eddy viscosity gradient $\nu_{T}^{\prime}=d\nu_{T}/dz$. In Sec. VI.2,
the wall-normal profiles of these two transfers are examined in order to
determine how they affect the eddy mode shapes. It is demonstrated in Sec.
VI.3 that artificially adjusting the eddy viscosity gradient can significantly
manipulate the eddy mode shapes, particularly in the near-wall region. Section
VI.4 discusses the types of interactions that are modelled by the eddy
viscosity.

###  VI.1 Eddy Dissipation

As explained in Ref. Symon  _et al._ (2021), the addition of eddy viscosity
introduces new dissipation terms into the kinetic energy balance for each
scale. These terms can be derived by expanding the viscous term in Eq. (9) and
Fourier-transforming in the homogeneous directions

$\bm{\nabla}\cdot\left[\nu_{T}(\bm{\nabla}\hat{\bm{u}}+\bm{\nabla}\hat{\bm{u}}^{T})\right]=\nu_{T}\bm{\nabla}^{2}\hat{\bm{u}}+(\bm{\nabla}\cdot\nu_{T})(\bm{\nabla}\hat{\bm{u}}+\bm{\nabla}\hat{\bm{u}}^{T}).$
(19)

Rewriting Eq. (19) in index notation, taking the inner product with respect to
$u_{i}$, and averaging the final expression in time yields

$\hat{\mathcal{D}}_{T}(k_{x},k_{y},z)=\underbrace{-\nu\overline{\frac{\partial\hat{u}_{i}}{\partial
x_{j}}\frac{\partial\hat{u}_{i}}{\partial
x_{j}}}}_{\hat{D}(k_{x},k_{y},z)}\underbrace{-(\nu_{T}(z)-\nu)\overline{\frac{\partial\hat{u}_{i}}{\partial
x_{j}}\frac{\partial\hat{u}_{i}}{\partial
x_{j}}}}_{\hat{V}(k_{x},k_{y},z)}\underbrace{+\frac{d\nu_{T}(z)}{dz}\overline{\hat{u}_{i}\left(\frac{\partial\hat{u}_{i}}{\partial
z}+\frac{\partial\hat{w}}{\partial x_{i}}\right)}}_{\hat{G}(k_{x},k_{y},z)},$
(20)

where $i,j=1,2,3$ and $\hat{\mathcal{D}}_{T}$ is dissipation due to molecular
and eddy viscosity. $\hat{\mathcal{D}}_{T}$ can be split into three separate
terms: $\hat{D}$ is dissipation due to molecular viscosity, $\hat{V}$ is
additional dissipation introduced by a wall-varying effective Reynolds number
$\nu_{T}(z)-\nu$, and $\hat{G}$ is dissipation due to the wall-normal gradient
of the eddy viscosity profile. The combined effect of $\hat{V}+\hat{G}$ can be
referred to as eddy dissipation, i.e. $\widehat{Edd}=\hat{V}+\hat{G}$. It can
be noted that without eddy viscosity, $\widehat{Edd}=\hat{V}=\hat{G}=0$ .

###  VI.2 Positive and negative energy transfers

It is guaranteed that $\hat{D}$ and $\hat{V}$ are real and negative at all
wall-normal locations, but this is not the case for $\hat{G}$. In fact, it can
be shown that $\hat{G}$ is likely to be positive near the wall. If the final
term of Eq. (20) is expanded in full, then $\hat{G}$ becomes

$\hat{G}=\frac{d\nu_{T}}{dz}\overline{\hat{u}\frac{\partial\hat{u}}{\partial
z}}+\frac{d\nu_{T}}{dz}\overline{\hat{v}\frac{\partial\hat{v}}{\partial
z}}+2\frac{d\nu_{T}}{dz}\overline{\hat{w}\frac{\partial\hat{w}}{\partial
z}}+\frac{d\nu_{T}}{dz}\overline{\hat{u}\frac{\partial\hat{w}}{\partial
x}}+\frac{d\nu_{T}}{dz}\overline{\hat{v}\frac{\partial\hat{w}}{\partial y}}.$
(21)

Assuming that the streamwise velocity component is significantly stronger than
the spanwise and wall-normal components, i.e. $\hat{u}\gg\hat{v},\hat{w}$, and
that wall-normal gradients dominate over streamwise and spanwise gradients,
Eq. (21) can be approximated as

$\hat{G}\approx\frac{d\nu_{T}}{dz}\overline{\hat{u}\frac{\partial\hat{u}}{\partial
z}}.$ (22)

For the lower wall, $\hat{u}=0$ due to the no-slip condition and both
$\hat{u}$ and $\frac{\partial\hat{u}}{\partial z}$ must be either both
positive or negative moving away from the wall. In either case, the product of
$\hat{u}$ and $\frac{\partial\hat{u}}{\partial z}$ is positive. As shown in
Fig. 16, $\frac{d\nu_{T}}{dz}>0$ near the lower wall so the product of all
three terms that appear in Eq. 22 is positive in the near-wall region.
Although several approximations have been made to arrive at this result, it is
significant because $\hat{V}$ and $\hat{G}$ attempt to model nonlinear
transfer processes. If $\hat{G}$ can be positive at a given wall-normal
location, then it implies that eddy viscosity can model locally both positive
and negative energy transfer. The net energy transfer due to eddy viscosity
that is obtained by integrating $\widehat{Edd}$ over the wall-normal direction
tends to be negative Symon  _et al._ (2021).

Figure 16: Wall-normal profiles of the eddy viscosity gradient
$\nu_{T}^{\prime}$ for various $Re_{\tau}$

Before examining profiles of $\hat{V}(z)$ and $\hat{G}(z)$ for specific
scales, it is worth commenting further on the eddy viscosity gradient profiles
in Fig. 16. Similar to the $\nu_{T}$ profiles in Fig. 13, it can be observed
that $Re_{\tau}$ has a negligible impact on $\nu_{T}^{\prime}$ for most wall-
normal locations. The maxima and minima of profiles become more extreme and
closer to the wall as $Re_{\tau}$ increases. This localization suggests that
$\hat{G}(z)$ is likely to be more concentrated in the near-wall regions and
have a larger magnitude for higher Reynolds numbers.

Figure 17(a) compares
$\hat{Q}(z)=-\overline{\hat{u}_{i}\frac{\partial}{\partial
x_{j}}\widehat{u_{i}u_{j}}}$, the time-averaged nonlinear transfer from DNS,
for $(k_{x},k_{y})=(0,4)$ with $\widehat{Edd}(z)$ in Fig. 17(b). Since the DNS
results are averaged in time, they include contributions from all frequencies
(wave speeds) whereas the eddy predictions are for the most energetic
frequency only. This is a reasonable approximation since $\omega=0$ dominates
over all other frequencies (recall that wave speed is ill-defined for
$k_{x}=0$ modes). There is good agreement between $\hat{Q}(z)$ and
$\widehat{Edd}(z)$ at nearly all wall-normal locations other than the near-
wall region. Eddy dissipation correctly predicts positive energy transfer near
the wall but its magnitude is too large. In Fig. 17(c), $\widehat{Edd}(z)$ is
split into the contributions from $\hat{V}(z)$ and $\hat{G}(z)$. Fig. 17(c)
reinforces that $\hat{G}(z)$ is responsible for the positive energy transfer,
as predicted by Eq. (22), and that $\hat{V}(z)$ accounts for the bulk of the
negative energy transfer.

Figure 17: (a) Time-averaged nonlinear transfer from DNS for
$(k_{x},k_{y})=(0,4)$ compared to (b) eddy dissipation and (c) its components
for $\omega=0$.

Figure 18(a) compares $\hat{Q}(z)$ for $(k_{x},k_{y})=(4,30)$ with
$\widehat{Edd}(z)$ in Fig. 18(b). The eddy predictions are for the most
energetic wave speed $c^{+}=10$. Good agreement can be observed between
$\hat{Q}(z)$ and $\widehat{Edd}(z)$ although the eddy dissipation is most
negative at $z^{+}=20$ instead of $z^{+}=15$. Similar to the previous scale,
eddy dissipation overestimates the positive energy transfer in the near-wall
region which, as seen in Fig. 18(c), is driven by $\hat{G}(z)$. The
contribution from $\hat{V}(z)$, on the other hand, is negative at all wall-
normal locations. It can be concluded that the eddy viscosity gradient plays
an important role in modeling positive energy transfer processes. The
predicted $\widehat{Edd}(z)$, nonetheless, exceeds the true nonlinear transfer
in the near-wall region. The next section considers the effect of artificially
weakening the eddy viscosity gradient to analyze its impact on the resulting
mode shapes.

Figure 18: (a) Time-averaged nonlinear transfer from DNS for
$(k_{x},k_{y})=(4,30)$ compared to (b) eddy dissipation and (c) its components
for $c^{+}=10$.

###  VI.3 Artificially adjusting the eddy viscosity gradient

In this section, the strength of the eddy viscosity term is altered by
introducing a scaling factor $\mathcal{G}$ such that Eq. (19) becomes

$\bm{\nabla}\cdot\left[\nu_{T}(\bm{\nabla}\hat{\bm{u}}+\bm{\nabla}\hat{\bm{u}}^{T})\right]=\nu_{T}\bm{\nabla}^{2}\hat{\bm{u}}+\mathcal{G}(\bm{\nabla}\cdot\nu_{T})(\bm{\nabla}\hat{\bm{u}}+\bm{\nabla}\hat{\bm{u}}^{T}),$
(23)

where $\mathcal{G}\geq 0$ controls the strength of the eddy viscosity gradient
term. Figure 19 illustrates the impact of artificially adjusting $\mathcal{G}$
for Mode 0. It can be seen that as $\mathcal{G}$ increases, the peak energy of
the streamwise and spanwise velocity components shifts closer to the wall. The
shape of the wall-normal velocity component is less affected but its magnitude
decreases. It can be reasoned that increasing $\mathcal{G}$ results in greater
positive energy transfer in the near-wall region. The energy of the modes,
consequently, is redistributed towards the wall and to the wall-parallel
velocity components. For some scales, $\mathcal{G}$ is too strong, resulting
in too much energy in the streamwise and spanwise velocity components near the
wall as seen for Mode 2 in Sec. V.1. These biased mode shapes explain why
linear-based estimation techniques that use an eddy viscosity model, e.g.
Refs. Madhusudanan  _et al._ (2019); Gupta  _et al._ (2021), overpredict the
strength of fluctuations in the near-wall region when measurements are known
in the logarithmic region. The biased mode shapes are also consistent with
Ref. Amaral  _et al._ (2021) who found that the eddy viscosity model
underestimates fluctuations when using wall-based measurements.

Figure 19: Leading eddy modes for $\bm{k}=(0,4,0)$ for various strengths of
the eddy viscosity gradient $\mathcal{G}$.

Figure 20 presents the effect of $\mathcal{G}$ on the Mode 1 shapes. Unlike
Mode 0, $\mathcal{G}$ has almost no impact on the structures although the
spanwise component becomes slightly less energetic relative to the other
velocity components as $\mathcal{G}$ increases. The main difference occurs in
the streamwise velocity component very close to the wall. The inset of Fig. 19
shows that the mode shape becomes flatter or more blunt as $\mathcal{G}$
increases. Similar to Mode 0, the flatter profiles can be attributed to
$\hat{G}$ injecting energy in this region of the flow. This also explains why
the constant eddy viscosity model, which has a wall-normal gradient of zero
everywhere, was not able to reproduce the near-wall behavior in Fig. 14(c)
that appeared in Fig. 14(b).

Figure 20: Leading eddy modes for $\bm{k}=(1,2,18.5)$ for various strengths of
the eddy viscosity gradient $\mathcal{G}$.

###  VI.4 Discussion

The previous sections have analyzed the energy transfer processes that can be
modeled by the Cess eddy viscosity profile. This section aims to contextualize
these results with respect to recent low-order modeling efforts in the
literature. One key challenge is to identify the smallest subset of nonlinear
interactions that are needed to sustain a wall-bounded turbulent flow at high
Reynolds numbers. The generalized quasi-linear (GQL) approximation Marston
_et al._ (2016), in particular, is able to reproduce key statistical features
of wall-bounded turbulence Hernández  _et al._ (2021). The GQL approximation
decomposes the flow into a low-wavenumber group and a high-wavenumber group.
Nonlinear interactions involving the high-wavenumber group are removed. If the
low-wavenumber group is restricted to the mean flow, the quasi-linear
approximation (QLA) is recovered Farrell and Ioannou (2007); Marston  _et al._
(2008); Thomas  _et al._ (2014, 2015).

Reference Hernández  _et al._ (2021) has shown that the GQL approach retains
triadic interactions that are responsible for the scattering mechanism and
inverse energy transfer in the near-wall region. The former can be attributed
to low-high wavenumber interactions that feed into high wavenumbers while the
latter arise from high-high wavenumber interactions that feed into low
wavenumbers. It is posited that the negative energy transfer modeled by eddy
viscosity through $\hat{V}$ can be interpreted as a scattering mechanism that
removes energy from larger scales and redistributes it to smaller scales. The
positive energy transfer modeled by $\hat{G}$, on the other hand, reproduces
the inverse energy transfer in the near-wall region. The combined effect of
these two transfers, therefore, encapsulates the effect of small scales on the
large scales.

##  VII Conclusions

The predictions of resolvent analysis with and without eddy viscosity have
been evaluated for a friction Reynolds number of $Re_{\tau}=550$. The accuracy
of the predictions were assessed using scalar measures including low-rank maps
of the operator, projection coefficients, and the most amplified wave speed as
well as direct comparison of the structures. The addition of eddy viscosity
distorted the linear amplification mechanisms that are identified by resolvent
analysis. As such, the low-rank maps of the eddy operator did not align well
with the turbulent kinetic energy spectrum in comparison to the low-rank maps
of the standard resolvent operator. For eddy analysis, low-rank behavior was
identified for structures with spanwise wavelengths of $\lambda_{y}^{+}\approx
80$ and $\lambda_{y}/h=3.5$ independent of the wave speed under investigation.
This is consistent with previous transient growth and optimal harmonic forcing
studies del Álamo and Jiménez (2006); Hwang and Cossu (2010).

Despite the distorted low-rank maps, the predictions from eddy analysis were
generally in better agreement with DNS in comparison to resolvent analysis.
The highest projection coefficients were obtained for eddy modes with wave
speeds of $c^{+}=10$ and $c^{+}=18.9$. These correspond to structures
associated with the near-wall cycle and structures that are most energetic at
$z/h=\pm 0.5$. The most amplified wave speed identified by eddy analysis was
also found to match the most energetic wave speed in DNS for these two types
of structures. For most wavenumber pairs, however, both resolvent and eddy
analysis failed to predict the most energetic wave speed. Resolvent analysis
overestimated the correct wave speed while eddy analysis underestimated the
correct wave speed.

The importance of wave speed on the projection coefficients motivated the
selection of specific scales for comparing the SPOD modes to the resolvent and
eddy modes. Consistent with previous studies, resolvent modes tended to be
localized around the critical layer McKeon and Sharma (2010) and the
streamwise velocity component was too strong relative to the spanwise and
wall-normal velocity components. All three velocity components of the eddy
modes, meanwhile, matched their SPOD counterparts for the wave speeds of
$c^{+}=10$ and $c^{+}=18.9$. For other wave speeds, eddy analysis mostly
underestimated the most energetic wave speed and predicted structures that
were energetic closer to the wall than the SPOD modes. It was concluded that
for higher Reynolds numbers, SPOD and eddy modes will agree for $c^{+}=10$ and
$c^{+}=U_{CL}^{+}-2$ due to the invariance of the Cess profile in the near-
wall region up to $z^{+}=15$, when it is premultiplied by $Re_{\tau}$, and in
the outer region where it is maximum at $z/h=\pm 0.5$. Modes for
$c^{+}=U_{CL}^{+}-2$ could also be obtained with reasonable accuracy from eddy
analysis by setting constant $\nu_{T}=Re_{T}=12.5$.

Finally, there was an investigation into the energy processes that are modeled
by the Cess profile. Two transfers were identified, one of which is an
effective Reynolds number that varies spatially. Similar to dissipation, the
effective Reynolds number term is guaranteed to be real and negative. The
second term originates from the eddy viscosity gradient, which tends to be
positive in the near-wall region. The combined effect, called eddy
dissipation, models well the actual nonlinear transfer taking place in DNS
although the positive energy transfer near the wall is overestimated. This
explains the tendency for eddy modes to be too energetic near the wall as
observed for higher Reynolds numbers in Ref. Symon  _et al._ (2020) and in
linear-based estimation, e.g. Ref. Madhusudanan  _et al._ (2019).

In order to further improve predictions from resolvent analysis, the Cess
profile could be replaced with a scale-dependent eddy viscosity. Reference
Pickering  _et al._ (2021), for example, solved an inverse problem to find the
eddy viscosity profile that maximized the projection of the leading resolvent
mode onto the leading SPOD mode. The drawback of this approach is that it
requires sufficient data in order to compute the leading SPOD mode. As such,
other approaches have modified the Cess profile using scaling arguments Gupta
_et al._ (2021) or a stochastic approach to model background turbulence that
can inject or dissipate energy of coherent waves Tissot  _et al._ (2021).
Regardless of the approach selected, the Cess profile provides a good initial
condition for optimizing an eddy viscosity or modeling the effect of nonlinear
terms as done in Ref. Zare  _et al._ (2017).

## Acknowledgments

The authors would like to acknowledge the use of the IRIDIS High Performance
Computing Facility, and associated support services at the University of
Southampton.

##  Appendix A Linear operators

After elimination of the pressure, the linearized Navier-Stokes equations can
be rewritten for the wall-normal velocity $\hat{w}$ and wall-normal vorticity
$\hat{\eta}=ik_{y}\hat{u}-ik_{x}\hat{v}$. The matrices $\bm{A}$, $\bm{B}$, and
$\bm{C}$ that appear in (4) are

$\bm{A}=\bm{M}\left[\begin{array}[]{cc}\mathcal{L}_{OS}&0\\\
-ik_{y}U^{\prime}&\mathcal{L}_{SQ}\end{array}\right],$ (24a)
$\bm{B}=\bm{M}\left[\begin{array}[]{ccc}-ik_{x}\mathcal{D}&-ik_{y}\mathcal{D}&-k^{2}\\\
ik_{y}&-ik_{x}&0\end{array}\right],$ (24b)
$\bm{C}=\frac{1}{k^{2}}\left[\begin{array}[]{cc}ik_{x}\mathcal{D}&-ik_{y}\\\
ik_{y}\mathcal{D}&ik_{x}\\\ k^{2}&0\end{array}\right].$ (24c)

Both $\mathcal{D}$ and ′ represent differentiation in the wall-normal
direction and $k^{2}=k_{x}^{2}+k_{y}^{2}$. The mass matrix $\bm{M}$ is defined
as

$\bm{M}(k_{x},k_{y})=\left[\begin{array}[]{cc}\Delta^{-1}&0\\\
0&\bm{I}\end{array}\right],$ (25)

where $\Delta=\mathcal{D}^{2}-k^{2}$ and $\bm{I}$ is the identity matrix. The
Orr-Sommerfeld $\mathcal{L}_{OS}$ and Squire $\mathcal{L}_{SQ}$ operators are

$\mathcal{L}_{OS}=-ik_{x}U\Delta+ik_{x}U^{\prime\prime}+(1/Re_{\tau})\Delta^{2},$
(26a) $\mathcal{L}_{SQ}=-ik_{x}U+(1/Re_{\tau})\Delta.$ (26b)

With the addition of eddy viscosity, they become

$\mathcal{L}_{OS}=-ik_{x}U\Delta+ik_{x}U^{\prime\prime}+\nu_{T}\Delta^{2}+2\nu_{T}^{\prime}\mathcal{D}\Delta+\nu_{T}^{\prime\prime}(\mathcal{D}^{2}+k^{2}),$
(27a) $\mathcal{L}_{SQ}=-ik_{x}U+\nu_{T}\Delta+\nu_{T}^{\prime}\mathcal{D}.$
(27b)

## References

  * McKeon and Sharma (2010) B. J. McKeon and A. S. Sharma, “A critical layer framework for turbulent pipe flow,” J. Fluid Mech. 658, 336–382 (2010).
  * Sharma  _et al._ (2016) A. S. Sharma, R. Moarref, B. J. McKeon, J. S. Park, M. D. Graham, and A. P. Willis, “Low-dimensional representations of exact coherent states of the Navier-Stokes equations from the resolvent model of wall turbulence,” Phys. Rev. E 93, 021102 (2016).
  * Moarref  _et al._ (2013) R. Moarref, A. S. Sharma, J. A. Tropp, and B. J. McKeon, “Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels,” J. Fluid Mech. 734, 275–316 (2013).
  * Skouloudis and Hwang (2021) N. Skouloudis and Y. Hwang, “Scaling of turbulence intensities up to ${R}e_{\tau}=10^{6}$ with a resolvent-based quasilinear approximation,” Phys. Rev. Fluids 6, 034602 (2021).
  * Luhar  _et al._ (2014) M. Luhar, A. S. Sharma, and B. J. McKeon, “Opposition control within the resolvent analysis framework,” J. Fluid Mech. 749, 597–626 (2014).
  * Nakashima  _et al._ (2017) S. Nakashima, K. Fukagata, and M. Luhar, “Assessment of suboptimal control for turbulent skin friction reduction via resolvent analysis,” J. Fluid Mech. 828, 496–526 (2017).
  * Toedtli  _et al._ (2019) S. S. Toedtli, M. Luhar, and B. J. McKeon, “Predicting the response of turbulent channel flow to varying-phase opposition control: Resolvent analysis as a tool for flow control design,” Phys. Rev. Fluids 4, 073905 (2019).
  * Symon  _et al._ (2021) S. Symon, S. J. Illingworth, and I. Marusic, “Energy transfer in turbulent channel flows and implications for resolvent modelling,” J. Fluid Mech. 911, A3 (2021).
  * Jin  _et al._ (2021) B. Jin, S. Symon, and S. J. Illingworth, “Energy transfer mechanisms and resolvent analysis in the cylinder wake,” Phys. Rev. Fluids 6, 024702 (2021).
  * Reynolds and Hussain (1972) W. C. Reynolds and A. K. M. F. Hussain, “The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments,” J. Fluid Mech. 54, 263–288 (1972).
  * Pujals  _et al._ (2009) G. Pujals, M. García-Villalba, C. Cossu, and S. Depardon, “A note on optimal transient growth in turbulent channel flows,” Phys. Fluids 21, 015109 (2009).
  * Hwang and Cossu (2010) Y. Hwang and C. Cossu, “Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow,” J. Fluid Mech. 664, 51–73 (2010).
  * Morra  _et al._ (2019) P. Morra, O. Semeraro, D. S. Henningson, and C. Cossu, “On the relevance of Reynolds stresses in resolvent analysis of turbulent wall-bounded flows,” J. Fluid Mech. 867, 969–984 (2019).
  * Kuhn  _et al._ (2021) P. Kuhn, J. Soria, and K. Oberleithner, “Linear modelling of self-similar jet turbulence,” J. Fluid Mech. 919, A7 (2021).
  * Morra  _et al._ (2021) P. Morra, P. A. S. Nogueira, A. V. G. Cavalieri, and D. S. Henningson, “The colour of forcing statistics in resolvent analyses of turbulent channel flows,” J. Fluid Mech. 907, A24 (2021).
  * Pickering  _et al._ (2021) E. Pickering, G. Rigas, O. T. Schmidt, D. Sipp, and T. Colonius, “Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets,” J. Fluid Mech. 917, A29 (2021).
  * Lumley (1967) J. L. Lumley,  _Stochastic tools in turbulence_ (Academic Press, 1967).
  * Picard and Delville (2000) C. Picard and J. Delville, “Pressure velocity coupling in a subsonic round jet,” Int. J. Heat Fluid Flow 21, 359–364 (2000).
  * Towne  _et al._ (2018) A. Towne, O. T. Schmidt, and T. Colonius, “Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis,” J. Fluid Mech. 847, 821–867 (2018).
  * Nogueira  _et al._ (2021) P. A. S. Nogueira, P. Morra, E. Martini, A. V. G. Cavalieri, and D. S. Henningson, “Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow,” J. Fluid Mech. 908, A32 (2021).
  * Amaral  _et al._ (2021) F. R. Amaral, A. V. G. Cavalieri, E. Martini, P. Jordan, and A. Towne, “Resolvent-based estimation of turbulent channel flow using wall measurements,” J. Fluid Mech. 927, A17 (2021).
  * Symon  _et al._ (2020) S. Symon, S. J. Illingworth, and I. Marusic, “Large-scale structures predicted by linear models of wall-bounded turbulence,” J. Phys. Conf. Series 1522, 012006 (2020).
  * Abreu  _et al._ (2020) L. I. Abreu, A. V. G. Cavalieri, P. Schlatter, R. Vinuesa, and D. S. Henningson, “Spectral proper orthogonal decomposition and resolvent analysis of near-wall coherent structures in turbulent pipe flows,” J. Fluid Mech. 900, A11 (2020).
  * Ellingsen and Palm (1975) T. Ellingsen and E. Palm, “Stability of linear flow,” Phys. Fluids 18, 487–488 (1975).
  * Landahl (1980) M. T. Landahl, “A note on an algebraic instability of inviscid parallel shear flows,” J. Fluid Mech. 98, 243–251 (1980).
  * Tissot  _et al._ (2021) G. Tissot, A. V. G. Cavalieri, and E. Mémin, “Stochastic linear modes in a turbulent channel flow,” J. Fluid Mech. 912, A51 (2021).
  * Lee and Moser (2019) M. Lee and R. D. Moser, “Spectral analysis of the budget equation in turbulent channel flows at high Reynolds numbers,” J. Fluid Mech. 860, 886–938 (2019).
  * Kawata and Tsukahara (2021) T. Kawata and T. Tsukahara, “Scale interactions in turbulent plane Couette flows in minimal domains,” J. Fluid Mech. 911, A55 (2021).
  * Doohan  _et al._ (2021) P. Doohan, A. P. Willis, and Y. Hwang, “Minimal multi-scale dynamics of near-wall turbulence,” J. Fluid Mech. 913, A8 (2021).
  * Hernández  _et al._ (2021) C. G. Hernández, Q. Yang, and Y. Hwang, “Generalised quasilinear approximations of turbulent channel flow: Part 2. Spanwise scale interactions,” arXiv:2112.01972v1 (2021).
  * Jovanović and Bamieh (2005) M. R. Jovanović and B. Bamieh, “Componentwise energy amplification in channel flows,” J. Fluid Mech. 534, 145–183 (2005).
  * Barthel  _et al._ (2021) B. Barthel, X. Zhu, and B. J. McKeon, “Closing the loop: nonlinear Taylor vortex flow through the lens of resolvent analysis,” J. Fluid Mech. 924, A9 (2021).
  * Reynolds and Tiederman (1967) W. C. Reynolds and W. G. Tiederman, “Stability of turbulent channel flow, with application to Malkus’s theory,” J. Fluid Mech. 27, 253–272 (1967).
  * Cess (1958) R. D. Cess, “A survey of the literature on heat transfer in turbulent tube flow,” Tech. Rep. (1958).
  * del Álamo and Jiménez (2006) J. C. del Álamo and J. Jiménez, “Linear energy amplification in turbulent channels,” J. Fluid Mech. 640, 5–26 (2006).
  * Welch (1967) P. Welch, “The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms,” IEEE. Trans. Audio Electroacoust. 15, 70–73 (1967).
  * Gibson  _et al._ (2019) J. F. Gibson, F. Reetz, S. Azimi, A. Ferraro, T. Kreilos, H. Schrobsdorff, M. Farano, A. F. Yesil, S. S. Schütz, M. Culpo, and T. M. Schneider, “Channelflow 2.0,” manuscript in preparation (2019).
  * Lee and Moser (2015) M. Lee and R. D. Moser, “Direct numerical simulation of turbulent channel flow up to ${R}e_{\tau}\approx 5200$,” J. Fluid Mech. 774, 395–415 (2015).
  * del Álamo and Jiménez (2003) J. C. del Álamo and J. Jiménez, “Spectra of the very large anisotropic scales in turbulent channels,” Phys. Fluids 15, L41–L44 (2003).
  * del Álamo  _et al._ (2004) J. C. del Álamo, J. Jiménez, P. Zandonade, and R. D. Moser, “Scaling of the energy spectra of turbulent channels,” J. Fluid Mech. 500, 135–144 (2004).
  * Gupta  _et al._ (2021) V. Gupta, A. Madhusudanan, M. Wan, S. J. Illingworth, and M. P. Juniper, “Linear-model-based estimation in wall turbulence: improved stochastic forcing and eddy viscosity terms,” J. Fluid Mech. 925, A18 (2021).
  * Hwang (2015) Y. Hwang, “Statistical structure of self-sustaining attached eddies in turbulent channel flow,” J. Fluid Mech. 767, 254–289 (2015).
  * Trefethen  _et al._ (1993) L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll, “Hydrodynamic stability without eigenvalues,” Science 261, 578–584 (1993).
  * Schmid and Henningson (2001) P. J. Schmid and D. S. Henningson,  _Stability and transition in shear flows_ (Springer, 2001).
  * Rosenberg and McKeon (2019) K. Rosenberg and B. J. McKeon, “Efficient representation of exact coherent states of the Navier-Stokes equations using resolvent analysis,” Fluid Dyn. Res. 51, 011401 (2019).
  * Madhusudanan  _et al._ (2019) A. Madhusudanan, S. J. Illingworth, and I. Marusic, “Coherent large-scale structures from the linearized Navier-Stokes equations,” J. Fluid Mech. 873, 89–109 (2019).
  * Vadarevu  _et al._ (2019) S. B. Vadarevu, S. Symon, S. J. Illingworth, and I. Marusic, “Coherent structures in the linearized impulse response of turbulent channel flow,” J. Fluid Mech. 963, 1190–1203 (2019).
  * Kuhn  _et al._ (2022) P. Kuhn, J. S. Müller, S. Knechtel, J. Soria, and K. Oberleithner, “Influence of eddy viscosity on linear modeling of self-similar coherent structures in the jet far field,” AIAA Scitech 2022 Forum , 0460 (2022).
  * Hwang (2016) Y. Hwang, “Mesolayer of attached eddies in turbulent channel flow,” Phys. Rev. Fluids 1, 064401 (2016).
  * Marston  _et al._ (2016) J. B. Marston, G. Chini, and S. M. Tobias, “Generalized quasilinear approximation: application to zonal jets,” Phys. Rev. Lett. 116, 214501 (2016).
  * Farrell and Ioannou (2007) B. F. Farrell and P. J. Ioannou, “Structure and spacing of jets in barotropic turbulence,” J. Atmos. Sci. 64, 3652–3665 (2007).
  * Marston  _et al._ (2008) J. B. Marston, E. Conover, and T. Schneider, “Statistics of an unstable barotropic jet from a cumulant expansion,” J. Atmos. Sci. 65, 1955–1966 (2008).
  * Thomas  _et al._ (2014) V. L. Thomas, B. K. Lieu, M. R. Jovanović, B. F. Farrell, P. J. Ioannou, and D. F. Gayme, “Self-sustaining turbulence in a restricted nonlinear model of plane Couette flow,” Phys. Fluids 26, 105112 (2014).
  * Thomas  _et al._ (2015) V. L. Thomas, B. F. Farrell, P. J. Ioannou, and D. F. Gayme, “A minimal model of self-sustaining turbulence,” Phys. Fluids 27, 105104 (2015).
  * Zare  _et al._ (2017) A. Zare, M. R. Jovanović, and T. T. Georgiou, “Colour of turbulence,” J. Fluid Mech. 812, 636–680 (2017).