1. Introduction
By linearising the Navier–Stokes equations around the turbulent mean profile, an array of analytical techniques can be employed for the study and modelling of turbulent flows (for recent reviews, see McKeon Reference McKeon2017; Jovanović Reference Jovanović2021). The success of this approach is largely due to the importance and robustness of linear amplification mechanisms within fully nonlinear dynamics. Analogously to the classical approach for modelling the mean profile with an eddy viscosity closure, an appropriate closure is required in the fluctuating velocity equations to replace the nonlinear term. This nonlinear term mediates energy transfer and interactions across different scales in turbulent flows. Consequently, to mimic the role of the nonlinear term played at specific length scales, a scale-dependent model of nonlinearity should be included in the linearised equations of motion.
To prescribe such a model, the role of the nonlinearity in transferring energy across scales needs to be well understood. Direct numerical simulation (DNS) precisely resolves nonlinear interactions (e.g. Kim, Moin & Moser Reference Kim, Moin and Moser1987; Hoyas & Jiménez Reference Hoyas and Jiménez2006; Lee & Moser Reference Lee and Moser2015). However, this remains extremely costly, particularly at Reynolds numbers of practical and theoretical interest. Efforts have been made to find suitable approximations for the nonlinear term. One such approach is a quasi-linear approximation, where the interactions through the nonlinear term are manipulated (for a recent review, see Marston & Tobias Reference Marston and Tobias2023). The velocity field is typically decomposed into large- and small-scale states, the exact definitions of which depend on the specific modelling objectives. By neglecting the nonlinear self-interactions of the small-scale state, the equations for the small-scale state are effectively linearised around the large-scale one. In the original formulations, where the large-scale state is defined to be time-averaged velocity in channel flow, the closure of the small-scale component was based on a marginal stability criterion (Malkus Reference Malkus1956). However, the turbulent mean profile in such a flow is typically linearly stable (Butler & Farrell Reference Butler and Farrell1992; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009). Therefore, more sophisticated closures are necessary, leading to a variety of modelling approaches in its modern variants. Examples include stochastic structural stability theory (Farrell & Ioannou Reference Farrell and Ioannou2007, Reference Farrell and Ioannou2012), direct statistical simulation (Marston, Conover & Schneider Reference Marston, Conover and Schneider2008; Tobias & Marston Reference Tobias and Marston2013), restricted nonlinear models (Thomas et al. Reference Thomas, Lieu, Jovanović, Farrell, Ioannou and Gayme2014, Reference Thomas, Farrell, Ioannou and Gayme2015; Bretheim, Meneveau & Gayme Reference Bretheim, Meneveau and Gayme2015; Farrell et al. Reference Farrell, Ioannou, Jiménez, Constantinou, Lozano-Durán and Nikolaidis2016), a quasi-linear approximation applied to exact coherent states (Pausch et al. Reference Pausch, Yang, Hwang and Eckhardt2019) and generalised quasi-linear approximations (Marston, Chini & Tobias Reference Marston, Chini and Tobias2016; Tobias & Marston Reference Tobias and Marston2017; Hernández, Yang & Hwang Reference Hernández, Yang and Hwang2021, Reference Hernández, Yang and Hwang2022; Luo, Hernández & Hwang Reference Luo, Hernández and Hwang2023).
While the Reynolds decomposition is used in a variety of situations, it may be employed for the quasi-linear framework, with the large-scale time-averaged mean state and the small-scale fluctuating velocity state. For wall-bounded turbulence, a quasi-linear approximation following the Reynolds decomposition was developed in Hwang & Eckhardt (Reference Hwang and Eckhardt2020), referred to as a minimal quasi-linear approximation (MQLA). In this MQLA, a self-consistent closure for the fluctuating velocity state was provided. Here and throughout the text, self-consistent refers to the fact that the Reynolds shear stress generated by the velocity fluctuations should be consistent with the Reynolds shear stress that appears in the mean momentum equation determining the mean profile. As the mean profile is empirically well approximated over a wide range of Reynolds numbers (e.g. Cess Reference Cess1958), the MQLA becomes a predictive framework over the range of Reynolds numbers. An extension has recently been provided, in a data-driven quasi-linear approximation (DQLA) (Holford, Lee & Hwang Reference Holford, Lee and Hwang2024), to account for the non-zero streamwise Fourier modes not included in the MQLA. To retain the predictive nature of the MQLA, the extension into the streamwise wavenumber domain was also made predictive. To this end, the self-similar nature of the energy-containing motions generated by the linearised Navier–Stokes equations and those present in real turbulent flow was exploited (Holford, Lee & Hwang Reference Holford, Lee and Hwang2023). A self-similar weight for the streamwise wavenumbers was developed for all Reynolds numbers, originally determined by matching the two-dimensional velocity spectra from the linearised Navier–Stokes equations to the DNS spectra at 
$Re_\tau \simeq 5200$ (
$Re_\tau$ is the friction Reynolds number). The resulting steady-state statistics and velocity spectra compare reasonably well with those of DNS, with Reynolds scaling behaviour consistent with available simulation and experimental data (Hwang & Eckhardt Reference Hwang and Eckhardt2020; Skouloudis & Hwang Reference Skouloudis and Hwang2021; Holford et al. Reference Holford, Lee and Hwang2024). The linearised fluctuation equations in the DQLA employ an eddy viscosity diffusion model and a forcing term to replace the nonlinearity. The role of this simple model has been described in numerous previous linear modelling studies (e.g. Reynolds & Hussain Reference Reynolds and Hussain1972; del Álamo & Jiménez Reference del Álamo and Jiménez2006; Hwang & Cossu Reference Hwang and Cossu2010; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019), but precisely how and why such a simple model effectively replaces the role of the nonlinearity involving highly complex interactions across a wide range of scales still needs to be answered. Recent studies by Symon, Illingworth & Marusic (Reference Symon, Illingworth and Marusic2021) and Symon et al. (Reference Symon, Madhusudanan, Illingworth and Marusic2023) provide insightful discussions on this issue by assessing the role of the eddy viscosity diffusion operator in resolvent analyses applied to the linearised fluctuation equations. In Symon et al. (Reference Symon, Illingworth and Marusic2021), the nonlinear transport in the minimal flow unit of near-wall turbulence, as well as in exact coherent states, was compared with that of the eddy viscosity diffusion model. An approximate budget was derived for the leading-order resolvent modes with integration over the wall-normal direction. It was demonstrated that without the eddy viscosity diffusion, the nonlinear term computed from DNS projects poorly over the leading resolvent modes of the linear model. This was reaffirmed in Symon et al. (Reference Symon, Madhusudanan, Illingworth and Marusic2023), where various predictions made with a resolvent analysis, with and without eddy viscosity, were compared with turbulent channel flow at 
$Re_\tau \approx 550$. They showed that the resolvent modes of the eddy viscosity enhanced resolvent operator provide a better basis, with findings similar to those of Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021). Furthermore, the eddy viscosity model provides mechanisms of wall-normal energy transfer (see also Hwang Reference Hwang2016), resulting in the resolvent mode structures comparing favourably with the spectral proper orthogonal decomposition (Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018) modes from DNS.
The aforementioned studies have provided valuable insights into the recent modelling efforts of turbulence using linearised fluctuation equations. However, these attempts are limited to low Reynolds numbers (e.g. Symon et al. Reference Symon, Illingworth and Marusic2021), at which the separation between inner and outer scales is very little, or are often restricted to a few wavenumber pairs and frequencies relevant to particular sets of energy-containing motions (e.g. Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021; Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023). The strength of the forcing inputs considered is often not self-consistent, hindering the understanding of the global features on the role of the eddy viscosity diffusion and the forcing in modelling of the nonlinear term, especially at high Reynolds numbers. In this respect, the self-consistent nature of the same form of the nonlinearity in the DQLA offers a way to assess the modelled nonlinear term at high Reynolds numbers. In particular, it will enable us to formulate a full spectral energy budget of the fluctuations generated by the linearised equations, which can be directly compared with that of DNS. With the forcing self-consistently determined in the DQLA, the role of the forcing and eddy viscosity diffusion operator can be clearly separated. This can then demonstrate the role of the eddy viscosity diffusion model in mimicking the nonlinear term at high Reynolds numbers.
The main objective of the present study is to analyse the energy transport mechanisms within the DQLA through a spectral budget analysis, which will be compared with a DNS counterpart from Lee & Moser (Reference Lee and Moser2019). In § 2, the spectral energy budget is first derived for the stochastic linear modelling frameworks used in previous studies (e.g. Farrell & Ioannou Reference Farrell and Ioannou1993; Jovanović & Bamieh Reference Jovanović and Bamieh2005; Farrell & Ioannou Reference Farrell and Ioannou2007). The budget in terms of the Orr–Sommerfeld–Squire system is also provided, and the turbulent transport is related to the modified form of the Lyapunov equation introduced in some previous studies (Zare, Jovanović & Georgiou Reference Zare, Jovanović and Georgiou2017; Abootorabi & Zare Reference Abootorabi and Zare2023). In § 3, the DQLA is briefly recapped. Particular care is taken in prescribing the model for the nonlinear term, where the additional terms in the budget present due to the inclusion of an eddy viscosity are given and discussed. In § 4, the results of the spectral budget analysis for the DQLA are compared with the DNS termwise. In particular, it is seen that the eddy viscosity term in the DQLA models the energy removal by the nonlinear term at integral length scales. Furthermore, the pressure–strain statistics reveal that the linear model for the fluctuation equations in the DQLA does not contain the streak instability mechanism, which plays an important role in redistributing energy from the streamwise component to the other components. A discussion of the entire budget and overall performance of the DQLA and the eddy viscosity model is given in § 5. Finally, concluding remarks are provided in § 6.
2. Lyapunov-like equations and the spectral energy budget
To start, a spectral budget equation is considered in the context of the Lyapunov-like equations. The computation of the white-in-time noise response to the linearised Navier–Stokes equations has been a standard practice (e.g. Bamieh & Dahleh Reference Bamieh and Dahleh2001; Jovanović & Bamieh Reference Jovanović and Bamieh2005; Madhusudanan, Illingworth & Marusic Reference Madhusudanan, Illingworth and Marusic2019; Jovanović Reference Jovanović2021; Holford et al. Reference Holford, Lee and Hwang2024). Therefore, relating the Lyapunov-like equations to the spectral energy budget equation used in, for example, Lee & Moser (Reference Lee and Moser2019) will enable us to establish a framework that compares a model based on linearised Navier–Stokes equations for fluctuations, such as the DQLA in Holford et al. (Reference Holford, Lee and Hwang2024), with DNS. In particular, the formulation of the DQLA in Holford et al. (Reference Holford, Lee and Hwang2024) is directly related to the spectral energy budget analysis of DNS data in Lee & Moser (Reference Lee and Moser2019), since the covariance formulated in the DQLA is computed by solving the Lyapunov equations in Holford et al. (Reference Holford, Lee and Hwang2024). From this, it will be seen that the spectral budget equation is a Lyapunov-like equation to be satisfied for each wavenumber pair, and it becomes a statistical budget of the spectral velocity covariance for the DQLA. No modelling assumptions are made for the nonlinear term here (additional modelling used in the DQLA is detailed in § 3, where the corresponding budget terms are also defined). Further details of the derivation are included in Appendix A (see also Appendix B for details on the computation of the pressure–velocity and forcing–velocity statistics). The intra-scale budget is then put into a form that can be compared with the budget equation statistics typically stored in DNS.
2.1. Turbulent channel flow
Fully developed incompressible turbulent flow is considered between two infinitely long and wide plates, with 
$x$, 
$y$ and 
$z$ denoting the streamwise, wall-normal and spanwise directions, respectively. The two parallel plates are separated by a distance of 
$2h$, with the lower and upper walls located at 
$y=0$ and 
$y=2h$. The velocity vector is denoted by 
$\boldsymbol{u} = (u,v,w)$ with components along the streamwise, wall-normal and spanwise directions, respectively. The velocity field is decomposed into a time-average mean state and a component fluctuating relative to this state, denoted by 
$\boldsymbol{U} = (U(y), 0 ,0)$ and 
$\boldsymbol{u}' = (u',v',w')$: i.e. the Reynolds decomposition is employed. The states are coupled through the following equations:
$$\begin{gather} \nu\frac{\textrm{d}U}{\textrm{d}\kern0.05em y} - \overline{u'v'} = \frac{\tau_w}{\rho}\left(1-\frac{y}{h}\right), \end{gather}$$
$$\begin{gather}\frac{\partial \boldsymbol{u}'}{\partial t} + (\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u}' + (\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{U} ={-}\frac{1}{\rho}\boldsymbol{\nabla} p' + \nu \nabla^2\boldsymbol{u}' + {\boldsymbol{\mathcal{N}}}, \end{gather}$$where
\begin{equation} {\boldsymbol{\mathcal{N}}} ={-}\boldsymbol{\nabla} \boldsymbol{\cdot}(\boldsymbol{u}'\boldsymbol{u}' - \overline{\boldsymbol{u}' \boldsymbol{u}'}). \end{equation}
Here an overbar denotes the time average, 
$p'$ is the fluctuating pressure, 
$\nu$ is the kinematic viscosity, 
$\rho$ is the density and 
$\tau _w$ is the shear stress at the wall. Equation (2.1a) is the mean momentum equation, which retains the nonlinear Reynolds shear stress term feeding back from the fluctuating state. The evolution of the fluctuating state remains fully nonlinear and exact for the purpose of this budget formulation. Given the homogeneous nature of the flow in the wall-parallel directions, the following Fourier transform for the states is considered:
\begin{equation} \hat{\boldsymbol{u}}'(t,y;k_x,k_z) = \iint_{-\infty}^{\infty} \boldsymbol{u}'(t,x,y,z)\,{\rm e}^{{\rm i}(k_x x + k_z z)}\,\textrm{d}\kern0.7pt x\, \textrm{d}z. \end{equation}The evolution equation for the fluctuating state is now rewritten in Fourier space, with the operator and state dependence on the wavenumber pair and wall-normal location dropped from the notation:
\begin{equation} \frac{\partial \hat{\boldsymbol{u}}}{\partial t} = \boldsymbol{A}_{\hat{\boldsymbol{u}}}\hat{\boldsymbol{u}} - \boldsymbol{G} \hat{p} + \hat{\boldsymbol{\mathcal{N}}}, \end{equation}with linear operators
\begin{equation} \boldsymbol{A}_{\hat{\boldsymbol{u}}} =\begin{bmatrix} \mathcal{L}_{SQ} & -(\mathcal{D}U) & 0 \\ 0 & \mathcal{L}_{SQ} & 0 \\ 0 & 0 & \mathcal{L}_{SQ} \end{bmatrix}, \end{equation}where
\begin{equation} \mathcal{L}_{SQ} ={-}\textrm{i}k_xU + \nu\varDelta \end{equation}is the Squire operator and
\begin{equation} \boldsymbol{G} = (1/\rho)[\textrm{i}k_x \quad \mathcal{D} \quad \textrm{i}k_z]^{\rm T} \end{equation}
is the gradient in column matrix operator form. Here 
$\mathcal{D}$ denotes differentiation along the wall-normal direction and 
$\varDelta = \mathcal {D}^2-k^2$ is the Laplacian in Fourier space, where 
$k^2 = k_x^2+k_z^2$. The boundary conditions are the no-slip conditions on the velocity state, 
$\hat {\boldsymbol{u}} = {\boldsymbol{0}}$ at 
$y = 0$ and 
$y = 2h$.
2.2. Lyapunov-like equations
To establish a local energy budget of (2.3a) in Fourier space, an evolution of the equation for the second-order statistics is derived here in terms of the Lyapunov-like equation. Following the formulation in Balakrishnan (Reference Balakrishnan1981) (see Appendix A for details and also Hoepffner et al. Reference Hoepffner, Chevalier, Bewley and Henningson2005; Zare et al. Reference Zare, Jovanović and Georgiou2017; Jovanović Reference Jovanović2021), the evolution of the covariance operator between velocity states at instances 
$t_1$ and 
$t_2$ is given by
\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \langle R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}(t_1,t_2)\hat{\boldsymbol{m}},\hat{\boldsymbol{n}}\rangle = \mathbb{E}\left[\left\langle \left. \frac{\partial \hat{\boldsymbol{u}}}{\partial t}\right\rvert_{t_1}, \hat{\boldsymbol{m}}\right\rangle\langle \hat{\boldsymbol{u}}(t_2),\hat{\boldsymbol{n}}\rangle^* + \langle \hat{\boldsymbol{u}}(t_1), \hat{\boldsymbol{m}}\rangle\left\langle \left.\frac{\partial \hat{\boldsymbol{u}}}{\partial t}\right\rvert_{t_2},\hat{\boldsymbol{n}}\right\rangle^*\right], \end{equation}
where 
$\langle \cdot, \cdot \rangle$ is the standard inner product:
\begin{equation} \langle \hat{\boldsymbol{m}}, \hat{\boldsymbol{n}}\rangle = \int_{0}^{2h} \hat{\boldsymbol{n}}^H \hat{\boldsymbol{m}}\,\textrm{d}\kern0.05em y. \end{equation}
Here, 
$(\cdot )^H$ is the complex conjugate transpose and 
$\hat {\boldsymbol{m}}(y; k_x, k_z)$ and 
$\hat {\boldsymbol{n}}(y; k_x, k_z)$ are arbitrary three-dimensional vectors. For steady-state statistics, the time derivative of the covariance operator has to vanish in the limit of 
$t\rightarrow \infty$. Setting the time lag to zero, 
$t=t_1 = t_2$, and taking the limit to (2.4) yields
\begin{equation} \lim_{t\to\infty}\frac{\textrm{d}}{\textrm{d}t} \langle R_{\hat{\boldsymbol{u }}\hat{\boldsymbol{u }}}(t,t)\hat{\boldsymbol{m}},\hat{\boldsymbol{n}}\rangle = \lim_{t\to\infty}\mathbb{E}\left[\left\langle \frac{\partial \hat{\boldsymbol{u}}}{\partial {t}}, \hat{\boldsymbol{m}}\right\rangle\langle \hat{\boldsymbol{u}},\hat{\boldsymbol{n}}\rangle^* + \langle \hat{\boldsymbol{u}}, \hat{\boldsymbol{m}}\rangle\left\langle \frac{\partial \hat{\boldsymbol{u}}}{\partial t},\hat{\boldsymbol{n}}\right\rangle^*\right] = 0. \end{equation}
Replacing 
$\partial {\boldsymbol{u}}/\partial t$ with (2.3a) then gives
\begin{equation} \lim_{t\to\infty}\langle \mathcal{L}_{b}(t,t)\hat{\boldsymbol{m}},\hat{\boldsymbol{n}} \rangle = 0, \end{equation}
where 
$\mathcal {L}_{b}(t,t)$ is the linear operator containing the time-dependent budget terms obtained by rearranging (2.6). For details, see Appendix A. For arbitrary 
$\hat {\boldsymbol{m}}$ and 
$\hat {\boldsymbol{n}}$ this leads to
\begin{align} \mathcal{L}_{b}^{\infty} &\equiv \lim_{t\to\infty}\mathcal{L}_{b}(t,t) \nonumber\\ &= \boldsymbol{A}_{\hat{\boldsymbol{u}}} R_{\hat{\boldsymbol{u }}\hat{\boldsymbol{u }}}^\infty + R_{\hat{\boldsymbol{u }}\hat{\boldsymbol{u }}}^\infty\boldsymbol{A}_{\hat{\boldsymbol{u}}}^{{{\dagger}}} -( \boldsymbol{G R}_{\hat{\boldsymbol{u}}\hat{p}}^\infty + R_{\hat{p}\hat{\boldsymbol{u}}}^\infty\boldsymbol{G}^{{{\dagger}}}) + ( R_{\hat{\boldsymbol{u}}\hat{\mathcal{N}}}^\infty + R_{\hat{\mathcal{N}}\hat{\boldsymbol{u}}}^\infty) = 0, \end{align}
where 
$(\cdot )^\infty$ denotes a steady-state operator.
For modelling purposes, this budget is typically formulated for the wall-normal velocity-vorticity state, 
$\boldsymbol{q} = [v'\ \eta ']^{\rm T}$, where 
$\eta ' = {\partial u'}/{\partial z} - {\partial w'}/{\partial x}$. In this formulation, the pressure is eliminated, and the resulting budget equation for the 
$\hat{\boldsymbol{q}}$ state is given by
\begin{equation} \mathcal{L}_{b,\hat{\boldsymbol{q}}}^{\infty} = \boldsymbol{A}_{\hat{\boldsymbol{q}}} R_{\hat{\boldsymbol{q}}\hat{\boldsymbol{q}}}^{\infty} + R_{\hat{\boldsymbol{q}}\hat{\boldsymbol{q}}}^{\infty} \boldsymbol{A}_{\hat{\boldsymbol{q}}}^{{{\dagger}}} + \boldsymbol{T}_{\hat{\boldsymbol{q}}} = 0, \end{equation}with
\begin{equation} \boldsymbol{T}_{\hat{\boldsymbol{q}}}=\boldsymbol{BR}_{\hat{\boldsymbol{q}}\hat{\boldsymbol{\mathcal{N}}}}^\infty + R_{\hat{\boldsymbol{\mathcal{N}}}\hat{\boldsymbol{q}}}^\infty \boldsymbol{B}^{{{\dagger}}}, \end{equation}
where 
$\mathcal {L}_{b,\hat{\boldsymbol{q}}}^{\infty } = \boldsymbol{B} \mathcal {L}_{b}^{\infty } \boldsymbol{B}^{{\dagger} }$ is the budget operator for 
$\hat{\boldsymbol{q}}$, 
$\boldsymbol{A}_{\hat{\boldsymbol{q}}}$ the Orr–Sommerfeld–Squire system and 
$\boldsymbol{B}$ a matrix operator which contains the definition of vorticity and a set of operations to remove the irrotational pressure field (see Appendices A and B). This is a Lyapunov-like equation, where 
$\boldsymbol{T}_{\hat{\boldsymbol{q}}}$ is not necessarily sign definite. In this case, the solution 
$R_{\hat{\boldsymbol{q}}\hat{\boldsymbol{q}}}^{\infty }$ is not guaranteed to be positive definite, as required. Such a modelling framework which relaxes the sign-definiteness 
$\boldsymbol{T}_{\hat{\boldsymbol{q}}}$ has recently been explored in Zare et al. (Reference Zare, Jovanović and Georgiou2017) and Abootorabi & Zare (Reference Abootorabi and Zare2023) (see also Georgiou (Reference Georgiou2002a) and Georgiou (Reference Georgiou2002b); Jovanović (Reference Jovanović2021) for a more rigorous derivation and review of the modelling framework, respectively). Alternatively, one can directly prescribe 
$\boldsymbol{T}_{\hat{\boldsymbol{q}}}$ to be positive definite and guarantee that 
$R_{\hat{\boldsymbol{q}}\hat{\boldsymbol{q}}}$ is positive definite through properties of the Lyapunov equation (Zhou, Doyle & Glover Reference Zhou, Doyle and Glover1996). Indeed, this is the case when the nonlinear term is replaced by a white-in-time forcing. The covariance of the model nonlinear term in this case is
\begin{equation}
R_{\hat{\boldsymbol{\mathcal{N}}}\hat{\boldsymbol{\mathcal{N}}}}(t_1,t_2)
= {R}_{\tilde{\boldsymbol{\mathcal{N}}}\tilde{\boldsymbol{\mathcal{N}}}}^{\infty}\delta(t_1-t_2),
\end{equation}
where 
$\delta (t_1-t_2)$ is the Dirac delta function and 
${R}_{\tilde {\boldsymbol{\mathcal {N}}}\tilde {\boldsymbol{\mathcal {N}}}}^{\infty }$ is the time-invariant operator giving the spatial covariance of 
$\mathcal {N}$. Using properties of the delta-correlated forcing, the model transport term becomes
\begin{equation} \boldsymbol{BR}_{\hat{\boldsymbol{q}}\hat{\boldsymbol{\mathcal{N}}}}^\infty + R_{\hat{\boldsymbol{\mathcal{N}}}\hat{\boldsymbol{q}}}^\infty \boldsymbol{B}^{{{\dagger}}} = \boldsymbol{BR}_{\tilde{\boldsymbol{\mathcal{N}}}\tilde{\boldsymbol{\mathcal{N}}}}^{\infty}\boldsymbol{B}^{{{\dagger}}}, \end{equation}giving the standard Lyapunov equation.
2.3. Spectral energy budget
The Lyapunov-like equations derived in § 2.2 essentially provide an energy budget equation for the velocity covariance operators to satisfy at a given wavenumber pair. As the covariance matrices are not typically stored in DNS, the statistics can be further reduced to the typical budget equation terms (e.g. Mansour, Kim & Moin Reference Mansour, Kim and Moin1988; Mizuno Reference Mizuno2016; Lee & Moser Reference Lee and Moser2019). In place of arbitrary functions in (2.7), 
$\hat {\boldsymbol{m}}$ and 
$\hat {\boldsymbol{n}}$, 
$\boldsymbol{e}_i\delta (y-y')$ and 
$\boldsymbol{e}_j\delta (y-y'')$ can be used, where 
$\boldsymbol{e}_i$ are the canonical basis vectors of 
$\mathbb {R}^3$. This provides the two-point wall-normal correlation statistics, 
$R_{\hat {\boldsymbol{u }}\hat {\boldsymbol{u }}}^\infty (y',y''; k_x, k_z)$. By setting 
$y'=y''=y$, a single wall-normal profile of statistics is then recovered, giving the typical energy budget statistics for DNS.
Below, the resulting energy budget equation is rearranged to be directly comparable with that in Lee & Moser (Reference Lee and Moser2019). The notation of the statistics follows Lee & Moser (Reference Lee and Moser2019), with 
${\mathsf{E}}_{ij}(y; k_x,k_z)$ denoting the one-point statistics of the Reynolds shear stress tensor at a wavenumber pair, where the subscript indices correspond to the components. The formulation of the spectral budget analysis here directly follows Lee & Moser (Reference Lee and Moser2019), to which the reader is referred for a more complete discussion, with the definitions consistent with other typical budget analyses (e.g. Mansour et al. Reference Mansour, Kim and Moin1988; Mizuno Reference Mizuno2016; Cho, Hwang & Choi Reference Cho, Hwang and Choi2018). The formulated budget equation is written in the following form:
\begin{equation} {\mathsf{E}}_{ij}^{P} + {\mathsf{E}}_{ij}^{\nu} + {\mathsf{E}}_{ij}^{\varPi} + {\mathsf{E}}_{ij}^{T} = 0, \end{equation}
for production, viscous, pressure and nonlinear (or turbulent) transport terms, respectively. Considering that the time-averaged mean profile consists of 
$(U(y), 0 ,0)$, the production term becomes
\begin{equation} {\mathsf{E}}_{ij}^{P} ={-}(\mathbb{E}[\hat{v}\hat{u}_j^*]\tilde{\delta}_{1i} + \mathbb{E}[\hat{u}_i\hat{v}^*]\tilde{\delta}_{1j})(\mathcal{D}U), \end{equation}
where 
$\tilde {\delta }_{ij}$ is the Kronecker delta. In the physical domain, this provides a source term for the turbulent kinetic energy in the streamwise component. The viscous term is written as
\begin{equation} {\mathsf{E}}_{ij}^{\nu} = \nu(\mathbb{E}[\hat{u}_j^*\varDelta \hat{u}_i] + \mathbb{E}[\hat{u}_i\varDelta \hat{u}_j^*]), \end{equation}which can then be split into a wall-normal transport term and viscous dissipation term:
\begin{equation} {\mathsf{E}}_{ij}^{\nu} = {\mathsf{E}}_{ij}^{\nu, D} + {\mathsf{E}}_{ij}^{\nu, \epsilon}, \end{equation}
where 
${\mathsf{E}}_{ij}^{\nu,D}$ and 
${\mathsf{E}}_{ij}^{\nu,\epsilon }$ satisfy
\begin{equation} \int_{0}^h {\mathsf{E}}_{ij}^{\nu,D}\,\textrm{d}\kern0.05em y = 0 \quad \forall \,k_x,\ k_z \end{equation}and
\begin{equation} {\mathsf{E}}_{ij}^{\nu,\epsilon} \leq 0 \quad \forall\, y, \ k_x,\ k_z. \end{equation}
In general, 
${\mathsf{E}}_{ij}^{\nu,D}$ has no net effect and simply transports turbulent kinetic energy across the wall-normal domain, with local positive source regions at a wall-normal location balanced by negative sink regions to satisfy this integral condition. The term 
${\mathsf{E}}_{ij}^{\nu,\epsilon }$ acts as a pure sink term, dissipating and removing turbulent kinetic energy, with
\begin{equation} {\mathsf{E}}_{ij}^{\nu,D} = \nu\mathcal{D}^2\mathbb{E}[\hat{u}_i\hat{u}_j^*] \end{equation}and
\begin{equation} {\mathsf{E}}_{ij}^{\nu,\epsilon} ={-}2\nu(k^2\mathbb{E}[\hat{u}_i\hat{u}_j^*] + \mathbb{E}[(\mathcal{D}\hat{u}_i)(\mathcal{D}\hat{u}_j)^*]). \end{equation}The correlation between the pressure gradient and velocity can be decomposed into a transport term and a traceless term:
\begin{equation} {\mathsf{E}}_{ij}^{\varPi} = {\mathsf{E}}_{ij}^{\varPi^s} + {\mathsf{E}}_{ij}^{\varPi^d}, \end{equation}
where 
${\mathsf{E}}_{ij}^{\varPi ^s}$ and 
${\mathsf{E}}_{ij}^{\varPi ^d}$ are
\begin{equation} {\mathsf{E}}_{ii}^{\varPi^s} = 0 \quad \forall \, y, \ k_x, \ k_z \end{equation}and
\begin{equation} \int_{0}^h {\mathsf{E}}_{ij}^{\varPi,d}\,\textrm{d}\kern0.05em y = 0 \quad \forall \,k_x,\ k_z. \end{equation}
There is no unique way of employing such a decomposition (Lumley Reference Lumley1975), and here the typical pressure–strain term is used for 
${\mathsf{E}}_{ij}^{\varPi,s}$, giving
\begin{equation} {\mathsf{E}}_{ij}^{\varPi, s} ={-}(\mathbb{E}[\hat{p}(\boldsymbol{G}_i \hat{u}_j)^*] + \mathbb{E}[\hat{p}^*\boldsymbol{G}_j \hat{u}_i] ), \end{equation}
where 
$\boldsymbol{G}_i$ is the 
$i$th entry in the gradient operator (2.3d), and the transport term becomes
\begin{equation} {\mathsf{E}}_{ij}^{\varPi, d} = \mathcal{D}\mathbb{E}[\hat{p}\hat{u}_j^*]\tilde{\delta}_{i2} + \mathcal{D}\mathbb{E}[\hat{p}^*\hat{u}_i]\tilde{\delta}_{j2}, \end{equation}with all transport occurring through the wall-normal velocity component. The traceless term by construction has zero net effect on the turbulent kinetic energy and does not appear in the budget for the total turbulent kinetic energy. However, for componentwise considerations, it acts as a source/sink term, dependent on the sign, with the other velocity components balancing this and satisfying the traceless condition.
Finally, the correlation between the nonlinear term and the velocity state, i.e. the turbulent transport term, can be decomposed. The total term is given by
\begin{equation} {\mathsf{E}}_{ij}^{T} = \mathbb{E}[\boldsymbol{G}_k\widehat{u_ku_i}(\hat{u}_j)^*] + \mathbb{E}[\hat{u}_i(\boldsymbol{G}_k\widehat{u_ku_j})^*]. \end{equation}In the spectral budget analysis of the DNS (Lee & Moser Reference Lee and Moser2019), this is split into an interscale transport term and a wall-normal transport term as follows:
\begin{equation} {\mathsf{E}}_{ij}^{T} = {\mathsf{E}}_{ij}^{T^\perp} + {\mathsf{E}}_{ij}^{T\parallel}, \end{equation}
where 
${\mathsf{E}}_{ij}^{T^\perp }$ satisfies the wall-normal transport condition, as in the previous decompositions, with
\begin{equation} \int_{0}^h {\mathsf{E}}_{ij}^{T^\perp}\,\textrm{d}\kern0.05em y = 0 \quad \forall \,k_x,\ k_z, \end{equation}which is given by
\begin{equation} {\mathsf{E}}_{ij}^{T^\perp} ={-}\tfrac{1}{2}\mathcal{D}(\mathbb{E}[\hat{u}_j(\widehat{u_iv})^*] + \mathbb{E}[\hat{u}_i^*(\widehat{u_jv})]). \end{equation}
The remaining term 
${\mathsf{E}}_{ij}^{T\parallel }$ (see Lee & Moser (Reference Lee and Moser2019) for the full expression) satisfies a condition that also ensures that 
${\mathsf{E}}_{ij}^{T\parallel }$ has zero net contribution to the turbulent kinetic energy. However, this is a global condition with
\begin{equation} \iint_{-\infty}^{\infty} {\mathsf{E}}_{ij}^{T^\parallel}\,\textrm{d}k_x \,\textrm{d}k_z = 0 \quad \forall \,y. \end{equation}
Hence, 
${\mathsf{E}}_{ij}^{T^\parallel }$ at a local wavenumber pair represents interscale transfer in the sense that turbulent transport is balanced out by the Fourier modes at other wavenumbers to satisfy this integral. That being said, both the wall-normal and interscale turbulent transport terms are precisely interscale, as they mathematically arise from triadic interactions across wavenumbers. For the model used in the present study, the nonlinear term is replaced by a forcing term and eddy viscosity diffusion operator (see § 3.1 for further details), and this condition is not designed to be satisfied, as is evident in § 4. No attempt is made to satisfy this global condition, as it is beyond the scope of this study. Instead, a separate decomposition for the nonlinear model used is provided in § 3.1.
Finally, with this budget formed in spectral space, taking the inverse Fourier transform, a budget in physical space is also satisfied:
\begin{equation} {\mathsf{P}}_{ij}(y) + {\mathsf{D}}_{ij}(y) + \epsilon_{ij}(y) + \varPi^S_{ij}(y) + \varPi^{\mathsf{D}}_{ij}(y) + {\mathsf{T}}^{}_{ij}(y) = 0, \end{equation}
with wall-normal profiles of production, viscous transport, dissipation, pressure–strain, pressure transport and turbulent transport, respectively. Note that due to (2.17d), the profile 
${\mathsf{T}}_{ij}(y)$ is determined solely from the 
${\mathsf{E}}_{ij}^{T^\perp }$ spectra in the DNS, while in the presented model, there are contributions from the non-conservative part of the nonlinear model.
3. Data-driven quasi-linear approximation
Having established the energy budget equation for the Navier–Stokes equations in Fourier space in terms of the Lyapunov-like equation in § 2, the DQLA is now introduced. Its modelling procedure is briefly discussed in § 3.1. The specific modelling feature of the DQLA (i.e. the nonlinear term model) is then related to the turbulent transport term (2.16) in the energy budget equation in § 3.2. Finally, the numerical solution procedure for the DQLA is briefly presented in § 3.3.
3.1. Formulation
This section provides a brief recap of the DQLA of Holford et al. (Reference Holford, Lee and Hwang2024) and the computational approach for the spectral energy budget formulated in the previous section. In the DQLA, the full nonlinear equation for the mean velocity is considered, while the equations for fluctuations are linearised about the mean. The model nonlinear term consists of an eddy viscosity diffusion term and a stochastic forcing term. In particular, for a given prescribed mean velocity, e.g. from Cess (Reference Cess1958), the stochastic forcing is determined such that the Reynolds shear stress from the fluctuating equations is identical to that from the mean equation, i.e. self-consistency. This feature allows us to directly compare the spectral energy budget of the DQLA with that of DNS, enabling us to fully understand the role of the eddy viscosity diffusion and stochastic forcing used in many previous linear modelling studies (e.g. Reynolds & Hussain Reference Reynolds and Hussain1972; Hwang & Cossu Reference Hwang and Cossu2010; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019; Symon et al. Reference Symon, Illingworth and Marusic2021).
The model for the nonlinear term in the DQLA is prescribed with the following form:
\begin{equation} \boldsymbol{\mathcal{N}} = \boldsymbol{\nabla} \boldsymbol{\cdot} (\nu_t(\boldsymbol{\nabla} \boldsymbol{u}' + (\boldsymbol{\nabla} \boldsymbol{u}')^T)) + \boldsymbol{f}, \end{equation}
with 
$\boldsymbol{f}$ being the stochastic forcing, and the eddy viscosity is frequently set to the empirical expression of Cess (Reference Cess1958):
\begin{equation} \nu_t =\frac{\nu}{2}\left[1+\frac{\kappa^2\,Re_\tau^2}{9}(1-y^2)^2(1+2y^2)^2 \left(1-\exp\left\{(|y|-1)\frac{Re_\tau}{A}\right\}\right)^2\right]^{{1}/{2}} - \frac{\nu}{2}. \end{equation}
In this study, the parameters were set to 
$A=28.7$ and 
$\kappa = 0.434$ based on a least-squares fit between the DNS mean profile at 
$Re_\tau \approx 5200$ and the mean profile determined with the eddy viscosity closure: 
$-\overline {u'v'} = \nu _t\,{\mathrm {d}U}/{\mathrm {d}\kern 0.05em y}$ in (2.1a). In Fourier space, the eddy viscosity diffusion term is given by
\begin{equation} \mathcal{L}_{\nu_t} = \begin{bmatrix} (\mathcal{D}\nu_t)\mathcal{D} + \nu_t\varDelta & {\rm i}k_x \mathcal{D}\nu_t & 0 \\ 0 & 2(\mathcal{D}\nu_t)\mathcal{D} + \nu_t\varDelta & 0 \\ 0 & {\rm i}k_z \mathcal{D}\nu_t & (\mathcal{D}\nu_t)\mathcal{D} + \nu_t\varDelta \end{bmatrix}. \end{equation}The stochastic forcing is designed to take the following form of the spectral covariance:
\begin{align} &\mathbb{E}[\hat{\boldsymbol{f}}(y,t;k_x,k_z)\hat{\boldsymbol{f}}^H(y',t;k_x,k_z) ] \nonumber\\ &\quad =\begin{bmatrix} W_u(k_x,k_z) & 0 & 0 \\ 0 & W_v(k_x,k_z) & 0\\ 0 & 0 & W_w(k_x,k_z) \end{bmatrix}\delta(y-y')\delta(t-t'), \end{align}
where 
$W_{r}$ are componentwise weights for the forcing with 
$r=\{u,v,w\}$.
Given the linear nature of the fluctuation equation for the DQLA, the spectral velocity covariance of the DQLA is written as
\begin{equation} \varPhi_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty} = \sum_{r=u,v,w}W_{r}(k_x,k_z)\varPhi_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}},r}^{\infty}, \end{equation}
where 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}}}^{\infty }=\mathbb {E}[\hat {\boldsymbol{u}}(y)\hat {\boldsymbol{u}}(y')^H]$ and 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty }$ is the spectral velocity covariance by solving the Lyapunov equation with componentenwise forcing: for example, 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},u}^{\infty }$ is obtained by setting 
$(W_u,W_v,W_w)=(1,0,0)$ (see Holford et al. (Reference Holford, Lee and Hwang2024) for further details). It was also shown that replacing 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty }$ with the one obtained from a few leading proper orthogonal decomposition (POD) modes improves the turbulence statistics from the DQLA, as this procedure effectively removes some of the unwanted non-physical features originating from 
$W_r$ not varying in the wall-normal direction (for a further discussion, see also Hwang & Eckhardt (Reference Hwang and Eckhardt2020) and Holford et al. (Reference Holford, Lee and Hwang2024)). For this purpose, the velocity covariance constructed with a few leading POD modes is further considered:
\begin{equation} \varPhi_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}},r}^{\infty,N_{POD}}(y,y';k_x,k_z) = \sum_{i}^{N_{POD}}\sigma_i \hat{\boldsymbol{u}}_{r,POD}(y,k_x,k_z)\hat{\boldsymbol{u}}_{r,POD}^H(y',k_x,k_z), \end{equation}
where 
$\hat {\boldsymbol{u}}_{r,POD}(y,k_x,k_z)$ is the POD mode obtained from 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty }$ and 
$\sigma _i$ the energy content of each POD mode with 
$\sigma _{i} \geq \sigma _{i+1}$. Consequently, the final form of the velocity covariance used in the DQLA is given as follows:
\begin{equation} \varPhi_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty} = \sum_{r=u,v,w}W_{r}(k_x,k_z)\varPhi_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}},r}^{\infty,N_{POD}}. \end{equation}
A value of 
$N_{POD} = 2$ is used, retaining the most energetic structure driven by the white-in-time forcing based on the symmetry in channel flow. This choice follows from Hwang & Eckhardt (Reference Hwang and Eckhardt2020), where the inclusion of higher-order POD modes introduces a non-physical peak in the velocity statistics towards the channel centreline due to the forcing in this region – this is also confirmed with the present DQLA by considering all the POD modes (not shown). For a discussion on the contributions of 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty }$ to 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}}}^{\infty }$, albeit with a white-in-time and wall-normal varying forcing rather than POD modes, the reader is referred to Holford et al. (Reference Holford, Lee and Hwang2023).
3.2. Spectral energy budget for DQLA
The only difference between the full Navier–Stokes equations and the equations for DQLA is the form of nonlinear term in (2.3a): i.e. (2.1c) versus (3.1). The difference in the spectral energy budget therefore appears only in the turbulent transport term. In the DQLA, this term is loosely referred to as a transport term, as it is a model for the turbulent transport term. However, strictly speaking, it is not a transport term, because the nonlinear term is not necessarily energy-preserving, given its form of an eddy-viscosity diffusion and forcing. Indeed, as is seen in § 4.6, its role only phenomenologically mimics that played by the nonlinear term in the full Navier–Stokes equations. In the DQLA, the turbulent transport term is written as
\begin{equation} {\mathsf{E}}_{ij}^{T} = {\mathsf{E}}_{ij}^{\nu_t} + {\mathsf{E}}_{ij}^{f}, \end{equation}
where 
${\mathsf{E}}_{ij}^{\nu _t}$ is the eddy-viscosity-related transport and
\begin{equation} {\mathsf{E}}_{ij}^{f}= \mathbb{E}[\hat{u}_j^*\hat{f}_i+\hat{f}_j^*\hat{u}_i] \end{equation}
is the forcing-related transport. The transport by the eddy viscosity operator can be further decomposed following 
${\mathsf{E}}_{ij}^{\nu }$ with
\begin{equation} {\mathsf{E}}_{ij}^{\nu_t} = {\mathsf{E}}_{ij}^{\nu_t, D} + {\mathsf{E}}_{ij}^{\nu_t, \epsilon} + {\mathsf{E}}_{ij}^{\nu_t, G}, \end{equation}
where 
${\mathsf{E}}_{ij}^{\nu _t,D}$ and 
${\mathsf{E}}_{ij}^{\nu _t,\epsilon }$ satisfy
\begin{equation} \int_{0}^{h} {\mathsf{E}}_{ij}^{\nu_t,D}\,\textrm{d}\kern0.05em y = 0 \quad \forall \,k_x,\ k_z \end{equation}and
\begin{equation} {\mathsf{E}}_{ij}^{\nu_t,\epsilon} \leq 0 \quad \forall\, y, \ k_x,\ k_z . \end{equation}
Hence overall the net effect of the eddy viscosity diffusion operator is the removal of energy through dissipation, given the positive profile of 
$\nu _t$, as discussed in Symon et al. (Reference Symon, Illingworth and Marusic2021, Reference Symon, Madhusudanan, Illingworth and Marusic2023). The 
${\mathsf{E}}_{ij}^{\nu _t, G}$ are the statistics that arise from the gradient of the wall-normal velocity component:
\begin{equation} {\mathsf{E}}_{ij}^{\nu_t, G} = (\mathcal{D}\nu_t)(\mathbb{E}[\hat{u}_i(\boldsymbol{G}_j v)^*] + \mathbb{E}[u_j^*\boldsymbol{G}_i \hat{v}]). \end{equation}These statistics are found to be small and largely negligible. The dissipation remains identical, with the eddy viscosity replacing the molecular viscosity:
\begin{equation} {\mathsf{E}}_{ij}^{\nu_t,\epsilon} ={-}2\nu_t(k^2\mathbb{E}[\hat{u}_i\hat{u}_j^*] + \mathbb{E}[(\mathcal{D}\hat{u}_i)(\mathcal{D}\hat{u}_j)^*]), \end{equation}
which remains a sink term, given the eddy viscosity profile remains positive across the wall-normal domain. The related transport terms vary due to the wall-normal variation of 
$\nu _t$ with
\begin{equation} {\mathsf{E}}_{ij}^{\nu_t,D} = {\mathsf{E}}_{ij}^{\nu_t,D_1} + {\mathsf{E}}_{ij}^{\nu_t,D_2}, \end{equation}where
\begin{equation} {\mathsf{E}}_{ij}^{\nu_t,D_1}=\nu_t \mathcal{D}^2\mathbb{E}[\hat{u}_i\hat{u}_j^*] \quad \textrm{and} \quad {\mathsf{E}}_{ij}^{\nu_t,D_2}=(\mathcal{D}\nu_t)\mathcal{D}\mathbb{E}[\hat{u}_i\hat{u}_j^*]. \end{equation}
These two terms are loosely referred to as transport terms in the current study, even though it is strictly their sum that satisfies (3.9b). They correspond to the wall-normal transport under 
$\nu _t$ and 
$\mathcal {D}\nu _t$, respectively, with similar terms offered in Symon et al. (Reference Symon, Madhusudanan, Illingworth and Marusic2023).
Finally, the inverse Fourier transform can be applied to the modelled turbulent transport term, yielding the turbulent transport profile in the wall-normal direction:
\begin{equation} {\mathsf{T}}_{ij}(y)={\mathsf{T}}_{ij}^{\nu_t}+{\mathsf{T}}_{ij}^f(y), \end{equation}where
\begin{equation} {\mathsf{T}}_{ij}^{\nu_t}(y) = {\mathsf{T}}_{ij}^{\nu_t,D_1}(y) + {\mathsf{T}}_{ij}^{\nu_t,D_2}(y) + {\mathsf{T}}_{ij}^{\nu_t,G}(y) + \epsilon_{ij}^{\nu_t}(y). \end{equation}
Here, 
${\mathsf{T}}_{ij}^{\nu _t,D_1}$, 
${\mathsf{T}}_{ij}^{\nu _t,D_2}$, 
${\mathsf{T}}_{ij}^{\nu _t,G}$ and 
$\epsilon _{ij}^{\nu _t}$ are from 
${\mathsf{E}}_{ij}^{\nu _t,D_1}$, 
${\mathsf{E}}_{ij}^{\nu _t,D_2}$, 
${\mathsf{E}}_{ij}^{\nu _t, G}$ and 
${\mathsf{E}}_{ij}^{\nu _t,\epsilon }$, respectively.
3.3. Numerical methods
The DQLA in this study was recently performed in Holford et al. (Reference Holford, Lee and Hwang2024). The wall-normal direction of the Orr–Sommerfeld–Squire system for each spatial wavenumber pair is discretised using a Chebyshev collocation method. The resulting discretised Lyapunov equation is solved using 
$\texttt {lyap}$ function in MATLAB to compute 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty }$. Then 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty,N_{POD}}(y,y';k_x,k_z)$ is subsequently obtained by computing appropriately weighted eigenvalues and eigenvectors of 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty }$ (Holford et al. Reference Holford, Lee and Hwang2024). With 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}},r}^{\infty,N_{POD}}(y,y';k_x,k_z)$ computed, the weight 
$W_{r}(k_x,k_z)$ is finally determined such that 
$\varPhi _{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}}}^{\infty }$ matches best with the two-dimensional velocity spectra from DNS at 
$Re_\tau \approx 5200$ (Lee & Moser Reference Lee and Moser2015), while the resulting Reynolds shear stress is numerically identical to that from the prescribed mean velocity given with (3.2). For this purpose, built upon the attached eddy hypothesis (Townsend Reference Townsend1976), the weight is decomposed into the following form:
\begin{equation} W_{r}(k_x,k_z)=W_{k_z}(k_z)W_{r,k_x}(k_x/k_z), \end{equation}
where 
$W_{k_z}(k_z)$ is a spanwise wavenumber dependent weight and 
$W_{r,k_x}(k_x/k_z)$ are self-similar streamwise weights for each of the velocity components. The weights are determined by solving the two related optimisation problems and for a more detailed account of the procedure, the reader is referred to Holford et al. (Reference Holford, Lee and Hwang2024). Since the DQLA formulation is based on the Orr–Sommerfeld–Squire system, the pressure–velocity and forcing–velocity correlations from the spectral energy budget analysis in §§ 2.2 and 3.2 need to be computed also. This is detailed in Appendix B. The Reynolds numbers for the DQLA considered are 
$Re_\tau = 1000$, 2000, 5200, 10 000, 20 000, and the numbers of grid points in the streamwise and spanwise directions for the different Reynolds numbers are fully documented in Holford et al. (Reference Holford, Lee and Hwang2024). Here, the number of wall-normal grid points was increased with 
$N_y = 256$, 386, 512, 768, 1024 for the Reynolds numbers considered in increasing order.
4. Results
In figure 1, the root-mean-square (r.m.s.) velocity profiles and the premultiplied budget terms between DNS and the DQLA are compared. A detailed discussion on the r.m.s. velocity profiles and their scaling behaviour with Reynolds number can be found in Holford et al. (Reference Holford, Lee and Hwang2024). Throughout, superscript ‘
$+$’ indicates the normalisation by the friction velocity 
$u_\tau$ and/or the viscous inner length scale 
$\delta _\nu = \nu /u_\tau$. Overall, the DQLA exhibits qualitative agreement with the DNS r.m.s. velocity statistics (figure 1a) across the components, the main missing feature being the plateau in the streamwise component. The scaling behaviour of the velocity profiles with Reynolds number is also consistent with DNS (omitted here; see Holford et al. (Reference Holford, Lee and Hwang2024) for further details), and is in agreement with the attached eddy hypothesis with finite Reynolds number corrections (Hwang, Hutchins & Marusic Reference Hwang, Hutchins and Marusic2022). The complete budget of the turbulent kinetic energy is presented in figure 1(b), taking the trace of (2.18), where the turbulent ‘transport’ of the DQLA consists of the sum of 
${\mathsf{T}}_{ii}^{\nu _t}$ and 
${\mathsf{T}}_{ii}^{f}$, as defined in (3.10). Up to 
$y^+ \approx 20$, the budgets exhibit qualitative agreement. Beyond this wall-normal location, only the dominant production term demonstrates similar behaviour, consistent with the DQLA model construction. In the DNS within this logarithmic region, production is effectively balanced by dissipation, with these two terms dominating budget considerations at a chosen wall-normal location. However, in the DQLA, production is approximately balanced by the model for the nonlinear term.
The wall-normal profiles of the turbulence statistics: (a) r.m.s. streamwise (black), wall-normal (blue) and spanwise (red) velocity fluctuations; (b) 
$y$-premultiplied turbulent kinetic energy budget consisting of production 
$P$ (black), dissipation 
$\epsilon$ (blue), viscous transport 
$D$ (red), pressure transport 
$\varPi ^D$ (green) and turbulent transport 
$T$ (yellow). Data obtained from DNS are the solid lines and data from the DQLA with 
$N_{POD}=2$ are the dashed lines. The logarithmic layer here appears approximately from 
$y^+=100$ to 
$y/h=0.15$ – for the mean velocity profiles of the DNS and DQLA, see figure 2 in Lee & Moser (Reference Lee and Moser2015) and figure 1 in Hwang & Eckhardt (Reference Hwang and Eckhardt2020), respectively.
4.1. Production
Figure 2 compares the premultiplied production profiles from DNS and the DQLA in inner- and outer-scaled wall-normal coordinates. There is a good agreement for the inner-scaled coordinates, with both DNS and the DQLA scaling well up to 
$y^+\approx 70$. The only small differences lie in the slightly lower peak in the DQLA, with this peak being broader: for example, 
$y^+P^+ \gtrapprox 5$ for 
$y^+ \in [10,60]$ as opposed to 
$y^+ \in [10,40]$ for the DNS. When comparing the profiles in outer-scaled coordinates, the DQLA scales much better, with well-defined outer scaling for 
$y/h \gtrapprox 0.1$. On the contrary, in the DNS a peak occurs in the logarithmic region, with a Reynolds-dependent magnitude and location, at least for the Reynolds numbers considered here. Despite this difference, the production of DQLA captures the important scaling behaviour for the logarithmic region, with a plateau in the production profile indicating inverse scaling with 
$y$. The magnitude of this plateau is consistent with the Cess profile of the eddy viscosity and its prescribed von Kármán constant, with the production in the logarithmic layer acting as
\begin{equation} y^+ P^+= {-}2 y^+ \overline{u'v'}^+\frac{\textrm{d}U^+}{\textrm{d}\kern0.05em y^+} \simeq 2y^+ \nu_t^+ \left(\frac{1}{\kappa^2y^+}\right)^2 = \frac{2}{\kappa}, \end{equation}
with 
$\nu _t \sim y$ being used, and 
${2}/{\kappa } \simeq 4.61$. Given this plateau occurs in the logarithmic region, the overall production there increases with 
$Re$ as the inner and outer length scales are better separated to define the logarithmic region. Overall, the comparison of the production profiles is related to the input of the DQLA model (i.e. the mean profile or the Cess eddy viscosity) rather than to the modelling framework itself, as this production profile is purely determined by the eddy viscosity closure used. Therefore, the framework could benefit from a more precise closure than the Cess profile used, particularly for modelling the outer peak and its scaling behaviour.
The 
$y$-premultiplied production profiles from (a,b) the DQLA and (c,d) DNS in (a,c) inner- and (b,d) outer-scaled wall-normal coordinates. Here, 
$Re_\tau \approx 1000$, 2000, 5200 for DNS and 
$Re_\tau = 1000$, 2000, 5200, 10 000, 20 000 for DQLA.
The premultiplied one-dimensional wavenumber spectra of production are compared in figure 3. The DQLA and DNS share the same qualitative features. There is a bimodal structure in the spanwise wavenumber spectra, with peaks occurring at approximately the same wall-normal locations and spanwise length scales in the near-wall and outer regions (i.e. the inner and outer peaks in figure 3b,d). Both the DQLA and DNS are energetic along a linear ridge, with the length scales exhibiting self-similar 
$y$ behaviour. In general, the DQLA is more energetic closer to the wall, with a length selection associated with long streaky motions, 
$y = 0.1 \lambda _z$ (figure 3b). On the other hand, the DNS is energetic closer towards the length scale associated with vortex packets, 
$y \approx 0.5 \lambda _z$ (Hwang Reference Hwang2015), with the spectra energetic between these two linear length scales (figure 3d). This trend is repeated in the streamwise wavenumber spectra with the DNS and DQLA (figure 3a,c). Again, the DQLA spectra are more energetic towards the associated streamwise length scale of the elongated streaky motions. It is also interesting to see that both DNS and the DQLA have a region of negative production in the streamwise spectra, with the peak occurring at approximately the same wall-normal location and streamwise length scale, although the DQLA region extends over the wall-normal domain and towards larger streamwise length scales. The origin of the negative production for the DQLA case is further explained in § 4.4. Here it is just pointed out that this region of negative production can be phenomenologically modelled by the forcing of the Orr–Sommerfeld–Squire system, considering that the DQLA spectra are constructed by weighting the most energetic structures of the linear model.
The 
$k_xy$- and 
$k_zy$-premultiplied one-dimensional production spectra from (a,b) DQLA and (c,d) DNS for (a,c) streamwise and (b,d) spanwise wavenumbers at 
$Re_\tau \approx 5200$. The solid lines are 
$y = 0.1 \lambda _z$ and 
$y = 0.01 \lambda _x$ and the dashed lines are 
$y = 0.5 \lambda _z$ and 
$y = 0.35 \lambda _x$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in 10 % increments.
4.2. Dissipation
Figure 4 compares the premultiplied dissipation profiles of the velocity components, with a lesser degree of agreement between DNS and the DQLA compared with the production profiles. The dissipation of the DQLA model is seen to present only in the near-wall region scaling in inner units, displaying very little values in the outer region (figure 4a–c). There is a well-defined near-wall peak at 
$y^+ \approx 8$ in the wall-parallel components, while the wall-normal component is considerably smoother, with a broad peak at 
$y^+ \approx 15$. The near-wall region in the wall-parallel components exhibits some Reynolds number dependence. In DNS (figure 4d–f), dissipation scales well in inner units up to 
$y^+ \approx 100$, although the amplitude for each component is seen to very mildly grow with 
$Re_\tau$. Below this location, anisotropic dissipation occurs, albeit with different characteristics when compared with the DQLA. In the DNS, the dissipation levels are roughly comparable in the streamwise perpendicular components, with the streamwise component dominating dissipation in the near-wall region. In the DQLA model, the two wall-parallel components more closely resemble each other, with the spanwise component having a reduced magnitude, while the wall-normal component has the least significance. The subtle Reynolds number dependence observed in the DNS spectra for the wall-parallel components in the near-wall region is notably weaker when compared with the DQLA model, with Reynolds dependence around the peak also. The primary qualitative differences between DNS and the DQLA lie in the absence of any features in the logarithmic/outer region in the DQLA model. In the DNS there is an approximate plateau with isotropic dissipation across the components, which provides an approximate balance for production in the logarithmic region. Given that the dissipation in the DQLA is present only in the near-wall region across different Reynolds numbers, another term in the budget must be balancing production, namely the nonlinear model in the DQLA from figure 1.
The 
$y$-premultiplied dissipation profile from (a–c) DQLA and (d–f) DNS in inner-scaled wall-normal coordinates for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components. Here 
$Re_\tau \approx 1000$, 2000, 5200 for DNS and 
$Re_\tau = 1000$, 2000, 5200, 10 000, 20 000 for the DQLA.
The premultiplied one-dimensional streamwise wavenumber spectra are compared in figure 5, with a similar comparison in the spanwise domain (not shown here). Overall, the spectra reflect the differences highlighted by the dissipation profiles: there is very weak/absent dissipation away from the wall, say for 
$y^+ \gtrapprox 100$ in all components for the DQLA (figure 5a–c). In particular, in the DQLA the dissipation in the logarithmic region occurs along a length scale much more local to production than in DNS, following a linear length scale selection slightly further from the wall than the linear length scale selection in production. In DNS the length scale selection for the strong dissipation region follows the local Kolmogorov length scale (figure 5d–f). In contrast, in the near-wall region, say for 
$y^+ \lessapprox 10$, the streamwise and spanwise dissipation spectra spread over a wide range of streamwise length scales for both DQLA and DNS, confirming the previously proposed mechanisms on their association with the wall-attached footprints observed in the streamwise and spanwise velocity spectra (Cho et al. Reference Cho, Hwang and Choi2018; Holford et al. Reference Holford, Lee and Hwang2024). However, some quantitative differences are also seen. The DQLA streamwise and spanwise dissipation spectra are relatively similar in their overall shape, with a near-wall peak and dissipation extending to the large-scale attached regions. Although the dissipation in the streamwise component of the DQLA may be reflective of the DNS for 
$y^+ \lessapprox 10$, the dissipation of the DQLA is much larger in the relative contribution to dissipation. Moreover, the dissipation of the spanwise component in the DNS is more similar to the wall-normal one, with the wall-attached features much more subtle.
The premultiplied streamwise wavenumber spectra for dissipation from (a–c) DQLA and (d–f) DNS for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components at 
$Re_\tau \approx 5200$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in increments of 10 %. The solid lines are 
$y/\lambda _x = 0.03$, 0.5 and 
$0.1$, and 
$y^{1/4}/\lambda _x = 13$, 20 and 
$20$ for streamwise, wall-normal and spanwise components in the DQLA and DNS, respectively.
4.3. Viscous transport
The viscous transport profiles are shown in inner-scaled wall-normal coordinates in figure 6, with a strong level of agreement between the DQLA and DNS. Both are highly anisotropic, with the streamwise component being the dominant term. The viscous transport shows weak Reynolds dependence, with this more pronounced in the DQLA. Three distinct near-wall peaks occur at 
$y^+ \approx 3,\ 13$ and 
$10$ that transfer energy towards and away from the wall, respectively. The most notable difference lies in the spanwise component, with this term in DNS being much less significant than the DQLA counterpart. The spanwise component shares the same qualitative features as the streamwise component at a reduced magnitude. This is true for both DNS and the DQLA, with the stronger viscous transport in the DQLA having a similar degree of Reynolds dependence to the streamwise component, with the Reynolds dependency in DNS being negligible relative to the streamwise component.
The 
$y$-premultiplied viscous transport profile from (a–c) the DQLA and (d–f) DNS in inner-scaled wall-normal coordinates for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components. Here 
$Re_\tau \approx 1000$, 2000, 5200 for DNS and 
$Re_\tau = 1000$, 2000, 5200, 10 000, 20 000 for the DQLA.
The premultiplied wavenumber spectra for the streamwise component of viscous transport are shown in figure 7, with strong agreement between DQLA and DNS. All three peaks are found at common streamwise and spanwise length scales, 
$(\lambda _x^+ ,\lambda _z^+) \approx (800,100)$ for DNS, and slightly larger scales in the DQLA, 
$(\lambda _x^+ ,\lambda _z^+) \approx (2000,200)$. Both the DNS and DQLA also show a region of viscous transport for the larger scales, accounting for the Reynolds dependency in the wall-normal profiles. This is further confirmed in the two-dimensional spectra, and is present in the DQLA also, with growth along the linear ridge for larger streamwise and spanwise length scales at this near-wall location (not shown). This Reynolds dependency in the near-wall region occurs largely for elongated streaky motions, similar to production, for say 
$\lambda _x/\lambda _z \gtrapprox 10$, with the large-scale motions becoming increasingly significant with Reynolds number.
The premultiplied spectra for viscous transport in the streamwise component from (a–c) the DQLA and (d–f) DNS for the (a,d) streamwise, (b,e) spanwise and (c, f) two-dimensional wavenumber spectra at 
$y^+ \approx 15$ for 
$Re_\tau \approx 5200$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in increments of 10 %. The solid lines are 
$\lambda _z/\lambda _x = 0.2$ and 
$0.1$ for the DQLA and DNS, respectively.
The strong similarity of viscous transport is due to the strong similarity in the velocity spectra between DQLA and DNS (Holford et al. Reference Holford, Lee and Hwang2024). Indeed, using the Taylor expansion, the near-wall two-dimensional viscous transport spectra may be written such that
\begin{equation} {\mathsf{E}}_{ij}^+=({\mathsf{E}}_{ij}^{\nu,D})^+{y^+}^2 +O({y^+}^4), \end{equation}
where the superscript 
$(\cdot )^+$ denotes the viscous inner scaling and 
${\mathsf{E}}_{ij}=\mathbb {E}[\hat {u}_i\hat {u}_j^*]$.
4.4. Pressure–strain
The premultiplied pressure–strain profiles are shown in figure 8. By construction, this term is traceless; therefore, the profiles can be interpreted as a statistical energy transfer between the components. Importantly, for components with significant source terms, for example, the production term in the streamwise component, the pressure–strain can redistribute the source term to other components. When comparing the profiles, the streamwise and spanwise share the same sign when comparing the DNS and DQLA, with the streamwise component negative (figure 8a,d), consistent with the redistribution of the production source term, and the spanwise component positive across the entire wall-normal domain (figure 8c, f). However, the streamwise component in the DQLA is at a significantly reduced magnitude when compared with the DNS, indicating that more of the production source term is retained in the streamwise component in the DQLA model. The wall-normal component in the DQLA is primarily negative, the opposite sign of the DNS (figure 8b,e). This can be associated with the fundamental limitations of the linear model used in the DQLA, where the streak instability/transient growth is not present by construction – these processes have been understood to redistribute the streamwise velocity fluctuations to cross-streamwise components in DNS (Cho et al. Reference Cho, Hwang and Choi2018; Doohan, Willis & Hwang Reference Doohan, Willis and Hwang2021). Here, in the DQLA model, the wall-normal velocity fluctuations arise primarily from the forcing, acting as a production-like term, explaining the differences in sign. It is also worth noting that this sign difference in the wall-normal component of pressure–strain is also present in the near-wall ‘splat’ region (Mansour et al. Reference Mansour, Kim and Moin1988; Lee & Moser Reference Lee and Moser2019) for 
$y^+ \lessapprox 10$, where the impermeability condition/continuity is enforced, shifting energy to the wall-parallel components. Moreover, in the DNS this ‘splat’ region is a transfer between the wall-normal and spanwise components, whereas in the DQLA the transfer is primarily between the streamwise and wall-normal fluctuations. Aside from the qualitative differences in the profiles, both DNS and DQLA appear to show a weak Reynolds number dependence in the near-wall region (streamwise and spanwise components, in particular). The profiles of the DQLA approximately collapse up to 
$y^+ \approx 12$, except in the streamwise component, compared with the DNS, which does so up to 
$y^+ \approx 50$. This difference between the DQLA and DNS is reflective of the differences in the production profiles (figure 2). The streamwise component of pressure–strain in the DQLA remains an approximate universal plateau, whereas in the DNS the peak in production in the logarithmic region (not modelled by the Cess profile) is presumably responsible for the similar behaviour in the DNS pressure–strain.
The 
$y$-premultiplied pressure–strain profile from (a–c) DQLA and (d–f) DNS in inner-scaled wall-normal coordinates for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components. Here 
$Re_\tau \approx 1000$, 2000, 5200 for DNS and 
$Re_\tau = 1000$, 2000, 5200, 10 000, 20 000 for the DQLA.
Figure 9 compares the premultiplied streamwise wavenumber pressure–strain spectra, which, like the profiles, are locally traceless for a considered wavenumber pair and sum to zero, and can be interpreted as a statistical intercomponent transfer. The spanwise wavenumber spectra compare similarly to the streamwise wavenumber spectra between the DNS and DQLA and are not shown, the only difference being in the streamwise component, which does not have a positive region in the spanwise wavenumber spectra. Both the DQLA and DNS exhibit pressure–strain spectra aligned along linear ridges (i.e. 
$\lambda _x \sim y$), with the length scales of the elongated streaky motion and vortex packets indicated by the solid and dashed lines, respectively, in the figures (Hwang Reference Hwang2015). In particular, for the DNS spectra (figure 9d–f), it becomes clear that the streak instability (or transient growth) mechanism is responsible for the intercomponent transfer with all the DNS spectra aligning well with 
$y \sim 0.35\lambda _x$ (Schoppa & Hussain Reference Schoppa and Hussain2002; de Giovanetti, Sung & Hwang Reference de Giovanetti, Sung and Hwang2017). There is little intercomponent transfer occurring for the elongated motions of 
$y \sim 0.01\lambda _x$, with production being retained in the streamwise component. In contrast, in the DQLA (figure 9a–c), there is no substantial mechanism for redistributing the turbulent kinetic energy from production across the components. In the streamwise pressure–strain component, the transfer is mainly towards the spanwise component, with the linear length scale aligning close to production, although slightly further from the wall (see figure 3). It is also worth pointing out that there is a region of positive pressure–strain in the streamwise component. This aligns well and is the same magnitude as the negative region in the streamwise wavenumber spectra of production (figure 3). This positive streamwise pressure–strain arises from the wall-normal component, indicating this negative production in the DQLA is a result of transfer from the wall-normal component to the streamwise one, with the wall-normal component itself arising from the forcing term, given its predominantly negative sign (see also Gupta et al. (Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021), where the modification of the forcing also removed the negative production). A similar feature is present in the spanwise component pressure–strain of DNS, with a negative region at relatively large streamwise wavelengths, with an intercomponent transfer from the wall-normal component. This negative region in the spanwise component and corresponding positive region in the wall-normal component are associated with streamwise vortices (Lee & Moser Reference Lee and Moser2019); with inspection of the two-dimensional spectra, this transfer is concentrated for streamwise-elongated structures. As expected, the physics of this is not present in the DQLA model, and this transfer is extremely weak upon inspection of the two-dimensional spectra (not shown here).
The 
$k_xy$-premultiplied streamwise wavenumber spectra for pressure–strain from (a–c) DQLA and (d–f) DNS for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components at 
$Re_\tau \approx 5200$. Contour levels are separated by 1 % of the maximum value of the streamwise component for DNS and spanwise component for the DQLA, increasing up to 10 % of the maximum, and then in 10 % increments. The solid lines are 
$y = 0.01\lambda _x$ and the dashed lines are 
$y = 0.35\lambda _x$.
4.5. Pressure transport
Next, the pressure transport term is compared. In Lee & Moser (Reference Lee and Moser2019), this term is grouped with the wall-normal turbulent transport, as both represent wall-normal transport arising from the nonlinearity. Here, the comparison is done separately, as the DQLA model has no constraints on pressure and the pressure–velocity correlation arising from the irrotational forcing is zero. Hence, grouping the pressure transport in the DQLA with the wall-normal nonlinear transport does not follow the same rationale.
The pressure transport is compared between the DQLA and DNS in figures 10 and 11. Here, there is little in common between the DNS and DQLA. The DNS pressure transport profile has a relatively simple structure, depositing energy in the near-wall region, with the source of this being further away from the wall (figure 10a). There is only weak Reynolds number dependence up to 
$y^+ \approx 30$, with approximately three inner scaling peaks at 
$y^+ \approx 7$, 
$20$ and 
$40$, and the profile eventually plateaus in the logarithmic region and becomes largely negative in the outer region. This is also reflected in the spanwise wavenumber spectra, where the spectra exhibit approximate linear scaling, removing energy away from the wall and depositing it towards the wall (figure 11a). The DQLA differs significantly from the pressure transport seen in DNS. For most of the wall-normal domain, the two have opposite signs, and the DQLA exhibits a much larger amplitude. In the transport profiles, there are still approximately three peaks at 
$y^+ \approx 7$, 
$20$ and 
$40$, however, with a sign change in pressure transport relative to the DNS (figure 10b), indicating that energy is mainly removed in the near-wall region and deposited in logarithmic/outer wall-normal locations. This is also reflected in the spectra, with most of the energy being transferred away from the wall along integral length scales (figure 11b). This is likely a consequence of neglecting any pressure-based constraints or irrotational forcing in the DQLA framework, with the focus solely on modelling velocity statistics.
The 
$y$-premultiplied pressure transport profiles from (a) DNS and (b) DQLA in inner-scaled wall-normal coordinates. Here 
$Re_\tau \approx 1000$, 2000, 5200 for DNS and 
$Re_\tau = 1000$, 2000, 5200, 10 000, 20 000 for the DQLA.
The 
$k_zy$-premultiplied pressure transport spanwise wavenumber spectra from (a) DNS and (b) DQLA at 
$Re_\tau \approx 5200$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in 10 % increments.
4.6. Turbulent transport and the nonlinear model
Before comparing the behaviour of the turbulent transport profile from the DNS and its model used in the DQLA, each term present in the DQLA nonlinear model is presented in figure 12. If the model (3.1) was energy conservative like the nonlinear term in the Navier–Stokes equations, the sum of all these terms integrated in the wall-normal domain should vanish. This condition is not imposed for the construction of the DQLA model (Holford et al. Reference Holford, Lee and Hwang2024), but the form of the model (3.1) allows it to have the mechanisms for both energy injection and removal through the forcing and eddy viscosity terms. Considering the wall-normal transport terms from the eddy viscosity integrated out to zero by the construction of the decomposition, the sum of the forcing ‘transport’ and dissipation is the injection/removal of energy by the nonlinear model at a given wall-normal location. Just from inspection of the profiles, it is clear that the model is non-conservative in each of the components. In the streamwise and spanwise components (figure 12a,c), the eddy viscosity dissipation results in the removal of energy for a majority of the wall-normal domain, with the only exception in the near-wall region of the spanwise component for 
$y^+ \lessapprox 10$ due to the forcing ‘transport’. In the wall-normal component (figure 12b), the forcing ‘transport’ is larger than the eddy viscosity dissipation, resulting in the injection of energy across the entire wall-normal domain. This is consistent with the observation in the pressure–strain spectra, where the wall-normal pressure–strain spectra are mainly negative, redistributing energy to the spanwise component.
The premultiplied nonlinear model profile from the DQLA decomposed into the 
$T^{f^s}$ (black), 
$-\epsilon ^{\nu _t}$ (blue), 
$T^{\nu _t,D_1}$ (red), 
$T^{\nu _t,D_2}$ (green) and 
$T^{\nu _t,G}$ (yellow) for the (a) streamwise, (b) wall-normal and 
$(c)$ spanwise components. In (b), 
$T_{22}^{\nu _t,D_2} = T_{22}^{\nu _t,G}$ and the plot is not repeated.
The eddy viscosity transport term, consisting of the ‘transport’ through eddy viscosity and the derivative of eddy viscosity, also represents the injection/removal of energy at a considered wall-normal location. As this term integrates out to zero across the wall-normal domain, the injection/removal at a given wall-normal location is balanced out by other regions, allowing the source location of local energy injections to be identified. Note that the ‘transport’ through the eddy viscosity profile (
$T^{\nu _t,D_1}$, the red line in figure 12) is the wall-normal viscous transport simply weighted by the eddy viscosity (see (3.9g)), with the same streamwise dominant term and with significant transfer in the spanwise component and negligible transfer in the wall-normal component. All three components of the transport by eddy viscosity are generally negative in the near-wall region while being positive away from the wall, indicating that this mechanism predominantly transfers the energy away from the wall. The transfer of energy in the proximity of the wall (
$y^+<5$) is generally highly damped compared with viscous transport (figure 6) by the low magnitude of the eddy viscosity. Moving a little away from the wall, the negative transfer of energy through the streamwise component becomes the main transport mechanism, with the peak at 
$y^+ \approx 11$. As the eddy viscosity profile substantially increases in the buffer layer (
$y^+ \approx 10\unicode{x2013}30$), the streamwise-component-dominant transport becomes increasingly substantial on increasing 
$y$. There is a peak in the positive transfer of energy for 
$y^+\approx 50$, with this now spanning the remainder of the wall-normal domain.
The ‘transport’ through the eddy viscosity derivative (
$T^{\nu _t,D_2}$, the green line) can similarly be regarded as the wall-normal derivative of the turbulence intensity profiles weighted by the wall-normal eddy viscosity gradient; see (3.9g). Taking into account that both the eddy viscosity and the velocity experience growth in the near-wall region due to the no-slip condition, this eddy viscosity derivative term is going to inject energy in the near-wall region (Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023). Moreover, the sign of the eddy viscosity derivative with the Cess profile remains positive for most of the wall-normal domain, as it monotonically increases with 
$\nu _t \sim y$, especially in the logarithmic region. Therefore, the overall sign of the transfer will depend on the derivative of the r.m.s. velocity profile. Given the impermeability condition at the wall, the wall-normal velocity spectra will always be negligible in this region, with the eddy viscosity derivative transfer term injecting energy primarily for the wall-parallel components. Consistent with the r.m.s. profiles, the eddy viscosity derivative ‘transport’ injects energy up to the near-wall peak in the r.m.s. velocity profiles at 
$y^+\approx 15$ in both the wall-parallel components. As 
$y$ increases, this transport becomes negative for both wall-parallel components, indicating that the injection of energy into the near wall by 
$T^{\nu _t,D_2}$ is a consequence of energy transfer towards the wall from the logarithmic and outer regions. This negative transport approximately balances the positive eddy viscosity transport 
$T^{\nu _t,D_1}$. Finally, it is noted that the statistics arising from the gradient of the wall-normal velocity component (the magenta lines in figure 12a,c and blue line in figure 12b) are largely negligible. The role of this term goes mainly into the pressure, particularly in the logarithmic region, where 
$(\mathcal {D}\nu _t) \sim 1$ and the gradient of the wall-normal velocity is irrotational, thereby not affecting the velocity field (see also the discussion in Appendix B).
With the role of the various terms in the nonlinear model established, the total eddy viscosity transport term (
$T^{\nu _t,D_1} + T^{\nu _t,D_2}$) is presented together, as well as the eddy viscosity dissipation and forcing transport term and compared with the turbulent transport from DNS in figure 13. The forcing transport and eddy viscosity dissipation profiles bear no resemblance to the turbulent transport in DNS. Instead, it becomes clear that the local balance of production and dissipation through molecular viscosity has been replaced by that of dissipation through eddy viscosity (see also figure 1b).
The premultiplied turbulent transport/nonlinear model profile from (a–f) DQLA and (g–i) DNS in inner-scaled wall-normal coordinates for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components. For DQLA, the transport profile is decomposed into (a–c) the forcing ‘transport’ and eddy viscosity dissipation and (d–f) the eddy viscosity transport. Here 
$Re_\tau \approx 1000$, 2000, 5200 for DNS and 
$Re_\tau = 1000$, 2000, 5200, 10 000, 20 000 for DQLA.
When comparing the scaling of the eddy viscosity transport profile (
$T^{\nu _t,D_1} + T^{\nu _t,D_2}$) in the DQLA (figure 13d–f) with the turbulent transport profile in DNS (figure 13g–i), it is remarkable to observe that there is strong agreement in both the qualitative features and the scaling behaviour, although no explanation can be given at this moment. In the streamwise component, the eddy viscosity transport well models the turbulent transport from DNS, with the positive and negative peaks at 
$y^+\approx 5$ and 
$y^+\approx 11$, respectively. Both DNS and the DQLA exhibit Reynolds-number-dependent peak magnitudes in these two near-wall transfers. The peak magnitudes are slightly increased in the DQLA, most noticeably the near-wall transfer towards the wall at 
$y^+\approx 5$. Moving away from the wall, the DQLA replicates the approximate positive logarithmic slope in the turbulent transfer from DNS for 
$y^+\gtrapprox 300$, although this emerges at higher Reynolds numbers in the DQLA (
$Re_\tau \gtrapprox 5200$) compared with the DNS (
$Re_\tau \gtrapprox 2000$). The other components compare less favourably. In the wall-normal component, the DQLA contains a negative peak in energy transfer at 
$y^+ \approx 100$, while in DNS this is much closer to the wall at 
$y^+ \approx 40$. Beyond this region, the DNS appears to scale well in inner units even up to the upper boundary of the logarithmic layer. There are even fewer common features in the spanwise component, with the DQLA transport being significantly dependent on Reynolds number.
The spanwise wavenumber spectra of the wall-normal turbulent transport from DNS and from the DQLA are shown in figure 14. In comparing them, it is clear that the eddy viscosity diffusion model well replicates the wall-normal turbulent transport from DNS, with the main transport being from the self-similar integral length scale (blue region), upwards away from the wall (dark red region). All three components also have weak positive transport in the region relatively close to the wall (the pale red region below blue region). However, in the wall-parallel components, these positive transport regions appear in the near-wall region. In the DQLA, this is relatively much stronger, particularly in the spanwise component. Moreover, in the spanwise component, this positive transport in the near-wall region follows the behaviour of the streamwise component, spanning outer- to inner-scaling spanwise length scales, while in the DNS, the transfer towards the wall mainly holds for the larger spanwise length scales.
The premultiplied spanwise wavenumber spectra for wall-normal turbulent transport from (a–c) DQLA and (d–f) DNS for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components at 
$Re_\tau \approx 5200$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in 10 % increments.
Finally, the spanwise wavenumber spectra of conservative interscale transport term from the DNS are compared with the non-conservative DQLA forcing and eddy viscosity dissipation terms in figure 15. As the DNS spectra are energy-conservative, the sink and source regions can be identified, and the related transport is across the different spanwise length scales. Here, there is the most significant difference between the DQLA and the DNS. In the DQLA (figure 15a–c), the nonlinear model mainly removes energy close to integral length scales, especially in the wall-parallel components, or acts as a driving term in the wall-normal case – the dissipation in the wall-normal component occurs slightly above this forcing ‘transport’ input, with the remnants of the eddy viscosity dissipation spectra present as the negative region in the spectra. As the eddy viscosity dissipation is mainly balancing production, the spectra also exhibit approximate 
$y$ scaling, indicating self-similarity. Although this self-similar scaling at the integral length scale in the logarithmic region is also present in the interscale transport of DNS (figure 15d–f), the only available mechanism for dissipation is through molecular viscosity. As such, the DNS has interscale transport through energy cascade down to the Kolmogorov length scale for dissipation, a feature not present in the DQLA. This is consistent with the observations on the dissipation in the DQLA, with the dissipation spectra following a linear length scale selection and being largely negligible in the logarithmic region. In both wall-parallel components, there is also a region of transport towards large scales in the near-wall region, originating from small scales (Cho et al. Reference Cho, Hwang and Choi2018; Lee & Moser Reference Lee and Moser2019). This is not present in the streamwise component of the DQLA. While these non-conservative spectra from the DQLA do not compare as favourably with the wall-normal transport by the eddy viscosity, the linear scaling feature in the eddy viscosity dissipation (and forcing in the wall-normal component) plays a crucial role in the energy balance (see the discussion in § 5.2).
The premultiplied spanwise wavenumber spectra for eddy viscosity dissipation and forcing transport summed together from the DQLA (a–c) and interscale turbulent transport from DNS (d–f) for the (a,d) streamwise, (b,e) wall-normal and (c, f) spanwise components at 
$Re_\tau \approx 5200$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in 10 % increments.
5. Discussion
The preceding section of results highlights the capabilities and limitations of the DQLA and its associated eddy-viscosity-enhanced linear model. The key points are summarised below:
(i) Production: When examining the sole source term in DNS, production, the DQLA demonstrates the most favourable comparison. The agreement in the production profile is expected due to the utilisation of the self-consistent Reynolds shear stress constraint in the DQLA. In particular, the DQLA accurately captures the plateau in the logarithmic region of the production profile, indicating the mean velocity proportional to
$y^{-1}$, consistent with the law of the wall. However, it misses the Reynolds-dependent peak observed in DNS, attributed to the simplistic Cess eddy viscosity profile employed. In terms of spectra, the self-similar behaviour of production spectra is evident, with streamwise and spanwise length scales following $y$ scaling, albeit closer to the wall than the DNS spectra. The success of this modelling is attributed to the application of the linearised Navier–Stokes equations along with the eddy viscosity model. Essentially, the DQLA enforces the attached eddy model, where the linearised Navier–Stokes equations, coupled with the eddy viscosity diffusion model, yield a self-similar structure in the wall-normal direction. These structures are superposed onto one another, guided by the self-consistent Reynolds shear stress constraint, resulting in self-similar premultiplied production spectra with most of the energetic content in streamwise-elongated structures.(ii) Componentwise energy redistribution: At the integral length scales, production is transferred across components via pressure–strain, with a negative streamwise pressure–strain that is comparable in magnitude to production. The DNS spectra suggest the streak instability as the mechanism for intercomponent transfer (Cho et al. Reference Cho, Hwang and Choi2018; Doohan et al. Reference Doohan, Willis and Hwang2021), demonstrated by well-aligned streamwise wavenumber spectra at
$y = 0.35\lambda _x$ (Hwang Reference Hwang2015). Through this mechanism, turbulent kinetic energy in the streamwise component primarily transfers to the wall-normal and spanwise components, serving as the primary statistical energy source at integral length scales in these components. In DNS, energy is also transferred from the spanwise to the wall-normal component through streamwise-elongated vortices, evident as the negative region in streamwise wavenumber spectra. Notably, this feature is absent in the DQLA. Production predominantly remains in the streamwise component, with the streamwise pressure–strain significantly lower than production. A wall-normal forcing term mimics the driving mechanism through the streak instability (or transient growth) and the subsequent nonlinear processes (Schoppa & Hussain Reference Schoppa and Hussain2002; de Giovanetti et al. Reference de Giovanetti, Sung and Hwang2017; Doohan et al. Reference Doohan, Willis and Hwang2021). Hence, this is only phenomenological, where the wall-normal component's pressure–strain is positive, transferring energy to the spanwise component. Despite differing signs in pressure–strain components, similar to production, the spectra display linear scaling, showcasing self-similar production transfer across components in DNS and redistributing the forcing in the DQLA.(iii) Energy transport and dissipation: Finally, with production redistributed across components, turbulent kinetic energy undergoes transport and dissipation. In DNS, the primary term that balances the remaining production in the streamwise component and pressure–strain in the spanwise and wall-normal components is the interscale turbulent transport term. This term removes energy at integral length scales approximately scaling with
$y$ in the logarithmic region, subsequently redistributing it towards the local Kolmogorov length scale through energy cascade under the energy-conservative nonlinear term in DNS. In contrast, the DQLA dissipates the source term locally in each component. The eddy viscosity dissipation serves as a substitute for interscale transport, removing energy at integral length scales and providing a self-similar energy balance, as well as replacing the role of molecular viscosity in balancing production globally in the wall-normal profiles. Despite neglecting energy cascade features, the self-consistent nature of the DQLA provides a sound energy balance, comparable to DNS and allows for a fair assessment of the eddy-viscosity-enhanced linear model. The eddy viscosity model contributes essential features to the linear model, which is discussed in the following subsections.
5.1. Eddy viscosity model as a colour of turbulence
Revisiting the problem formulation, one of the key features that the eddy viscosity diffusion provides is a colour of turbulence term in the same sense as in Zare et al. (Reference Zare, Jovanović and Georgiou2017). When examining the turbulent transport statistics from DNS, it becomes evident that a coloured-in-time model for the nonlinear term is indispensable for any linear model aiming to replicate the energy balance. In the white-in-time setting, where the nonlinear term is replaced solely by a forcing term delta-correlated in time, the turbulent transport term is equivalent to
\begin{equation}
R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{\mathcal{N}}}}^{\infty}+
R_{\hat{\boldsymbol{\mathcal{N}}}\hat{\boldsymbol{u}}}^{\infty} =
\tilde{R}_{\tilde{\boldsymbol{\mathcal{N}}}\tilde{\boldsymbol{\mathcal{N}}}}^{\infty},
\end{equation}
similar to (2.11), where 
$\tilde {\mathcal {N}}$ is a white-in-time forcing term and 
$\tilde {R}_{\tilde {\boldsymbol{\mathcal {N}}}\tilde {\boldsymbol{\mathcal {N}}}}^{\infty }$ contains the spatial statistics. Notably, as this term is an autocorrelation, the resulting turbulent transport statistics will strictly have a positive wall-normal profile. This poses important limitations since the primary role of interscale turbulent transport is to remove energy from integral length scales and redistribute it for dissipation. Therefore, in the white-in-time model, the statistics are inherently non-physical. With the introduction of a coloured-in-time model, the nonlinear term can now relax the sign definiteness issue for the white-in-time forcing. The resulting turbulence transport statistics are no longer restricted to be positive throughout the wall-normal domain. This highlights the essential nature of a coloured-in-time term in accurately replicating turbulent transport statistics, allowing for the removal and redistribution of energy across the wall-normal domain. Moreover, in the white-in-time setting, a conservative model for the nonlinear term is not possible from an interscale perspective. The statistics will always integrate out to be positive, given the model transport term is the variance profile of the white-noise forcing for each wavenumber.
While the preceding discussion advocates the use of a coloured-in-time model for the nonlinearity, the use of the eddy viscosity diffusion term in providing this can now be made more precise. Following the interpretation of this coloured-in-time term offered by Zare et al. (Reference Zare, Jovanović and Georgiou2017), the eddy viscosity diffusion term provides the colour of the turbulence model, that is, a modification to the linear dynamical operator. This modification preserves linear stability, allowing a white-in-time forcing to be used, where the Lyapunov equation of the modified linear operator can be solved to determine statistics. In the case of the eddy viscosity operator, following the Orr–Sommerfeld system, the modification of the linear dynamical operator gives the following equation for 
$\hat{\boldsymbol{q}}$:
\begin{equation} \frac{\partial \hat{\boldsymbol{q}}}{\partial t} = (\boldsymbol{A}_{\hat{\boldsymbol{q}}} + \boldsymbol{B} \boldsymbol{C}_f) \hat{\boldsymbol{q}} + \boldsymbol{B} \hat{\boldsymbol{f}}, \end{equation}where
\begin{equation} \boldsymbol{C}_f = \mathcal{L}_{\nu_t}\boldsymbol{C} \end{equation}
is the modification of the linear dynamical operator. Despite this relatively simplistic eddy viscosity term as a colour of turbulence model and it not offering the same efficacy in reproducing velocity statistics when compared with a more precise data-driven form (Zare et al. Reference Zare, Jovanović and Georgiou2017; Abootorabi & Zare Reference Abootorabi and Zare2023), it yields the advantage of having an explicit expression. This allows the eddy viscosity diffusion as the colour of turbulence model to be more easily interpreted. Furthermore, it retains a physical 
$\boldsymbol{B}$ operator in (A4), which effectively filters out the irrotational part of the nonlinear model (see also Appendix B), rather than a data-driven alternative in Zare et al. (Reference Zare, Jovanović and Georgiou2017). This colour of turbulence term is componentwise negative definite. It strictly removes energy through the eddy viscosity dissipation, with a conservative wall-normal transport term. With the DQLA framework self-consistently weighting the model to essentially match production, the statistics obtained from this simple yet effective colour of turbulence model can be compared with DNS, providing implications and potential improvements for the current eddy viscosity model.
As the problem formulation establishes a clearer connection between statistics from the Lyapunov-like equation for the linearised fluctuation model and the turbulent transport statistics from DNS, some immediate enhancements of the eddy viscosity model could be pursued by tuning it to match the DNS output – such a study was recently carried out for a turbulent jet by Pickering et al. (Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021) from a context of improving their resolvent analysis. Note that there is no reason for the eddy viscosity model in the fluctuating equation to be identical to that of the mean momentum equation. However, the parametrised form should probably possess the same key features as the Cess viscosity profile: (1) 
$\nu _t \sim y$ in the logarithmic region; (2) the same Reynolds-number scaling; (3) a wall-normal derivative to phenomenologically drive the footprints of the structures in the logarithmic and outer region (see figure 14).
First, at a given wavenumber pair, an eddy viscosity profile could be parametrised and tuned such that the wall-normal transport and eddy viscosity dissipation best match the associated features of the turbulent transport statistics. However, modifying the eddy viscosity profile will impact the DQLA framework, as the self-consistent forcing and weighting would need to change relative to the eddy viscosity profile, resulting in a coupled problem. To address this issue, the eddy viscosity could be optimised through a framework such as that of Holford et al. (Reference Holford, Lee and Hwang2023), where a white-in-time forcing of the eddy viscosity model is optimised to best match the velocity spectra from DNS. This could also allow a scale-dependent eddy viscosity profile to be used, given the qualitative nature of the interscale and wall-normal turbulent transport in DNS differs most from the current eddy viscosity model. Second, it is also clear that even though the application of forcing to the linearised Navier–Stokes equations yields anisotropic turbulence statistics due to the production term present only in the streamwise component, the anisotropy across the components does not match that of the DNS. Given that the main balance at integral scales essentially arises from the eddy viscosity term, the isotropic eddy viscosity diffusion term could also be tuned componentwise to reflect this. The dissipation, viscous transport and eddy viscosity wall-normal transport results suggest that the spanwise component would benefit most. In the DQLA, the response of the spanwise component is overly energetic and follows the same structure as the streamwise component. Whereas in the DNS the elongated streaky structures are relatively less significant in the spanwise component.
Lastly, the implications of using an eddy viscosity diffusion for colouring white noise (or even coloured) forcing and the possibility of a globally conservative model for the nonlinear term are discussed. The eddy viscosity model with the forcing should possess all the statistical features necessary for such a model: the eddy viscosity diffusion provides energy removal at integral scales and can mimic the turbulent wall-normal transport, while the white-in-time forcing could be tuned to deposit the energy at the Kolmogorov scales for dissipation under molecular viscosity. In fact, this is exactly the role of the energy cascade mediated by the nonlinear term in the original turbulent fluctuation equations. Therefore, the enforcement of global energy conservation is certainly possible using an eddy viscosity diffusion accompanied with a forcing from a technical viewpoint. It is, however, evident that it will not be able to replicate the dynamical processes involved in energy cascade, as the model is expected to remain only phenomenological. Such a modelling effort would not be largely beneficial to obtain novel physical insights. For the exploration of the dynamical processes of energy cascade and scale interactions, the generalised quasilinear approximations (e.g. Marston et al. Reference Marston, Chini and Tobias2016) would probably be more useful, as recently demonstrated by Hernández et al. (Reference Hernández, Yang and Hwang2021, Reference Hernández, Yang and Hwang2022) for turbulent channel flow.
5.2. Local self-similar energy balance
Despite the eddy viscosity diffusion model being limited to mimicking only wall-normal turbulent transport and energy removal at integral length scales, it will now be demonstrated that this eddy viscosity dissipation is essential for achieving a local self-similar budget. Following Hwang & Lee (Reference Hwang and Lee2020), this means the eddy viscosity diffusion model itself is consistent with a logarithmic mean profile, a generic feature that should be achieved by colour of turbulence models following these arguments. First, the mean-momentum equation in the logarithmic layer is revisited. At a sufficiently high Reynolds number, in the logarithmic region, the mean-momentum equation (2.1a) reduces to
\begin{equation} -\frac{\overline{u'v'}}{u_\tau^2} = 1. \end{equation}In terms of the two-dimensional spectra, this may be rewritten as
\begin{equation} \int_{-\infty}^{\infty} \int_{\infty}^{\infty}\frac{E_{12}}{u_\tau^2}\,\textrm{d}k_x\, \textrm{d}k_z ={-}1. \end{equation}
Given that the only relevant length scale in the logarithmic region is 
$y$, with 
$1/y\sim k_x \sim k_z$, self-similarity variables in the streamwise and spanwise directions are introduced through
\begin{equation} \xi_x=k_x y \quad \textrm{and} \quad \xi_z = k_z y. \end{equation}This results in (5.4) reducing to
\begin{equation} \int_{-\infty}^{\infty} \int_{\infty}^{\infty}\frac{E_{12}}{u_\tau^2y^2}\,\textrm{d}\xi_x \,\textrm{d}\xi_z ={-}1, \end{equation}which implies the following must be true:
\begin{equation} \frac{E_{12}}{y^2}={f(\xi_x,\xi_z)} \ \textrm{or, equivalently}\ {k_z^2E_{12}}=\frac{f(\xi_x,\xi_z)}{\xi_z^2}\equiv g(\xi_x,\xi_z). \end{equation}
The relation (5.7) results solely from the mean-momentum equation and the assumption that the only relevant length scale is 
$y$. Note that other self-similar variants for (5.7) are possible (e.g. using 
$k_x$) as long as the premultiplication factor for 
$E_{12}$ is dimensionally proportional to 
$1/y^2$ using the self-similar variables in (5.5). While (5.7) must be true in real flows, there is no guarantee the output of a linearised model for fluctuations satisfies this condition. This is precisely what the DQLA implements, albeit over the entire wall-normal domain, by ensuring the spectra are self-consistent with the mean shear present in the mean-momentum equation, with
\begin{equation} \int_{-\infty}^{\infty} \int_{\infty}^{\infty}E_{12}\,\textrm{d}k_x \,\textrm{d}k_z = \overline{u'v'}. \end{equation} With this theoretical property of the real flow spectra established and achieved in the linearised model by the DQLA self-consistent framework, now consider the fluctuating equations at a given spatial wavenumber pair. This can then be rewritten in terms of self-similarity variables. For convenience, take 
$k_z (\neq 0)$ as a given reference length scale (
$l = {1}/{k_z}$). First, the 
$\boldsymbol{A}_{\hat {\boldsymbol{u}}}$ operator is decomposed into the streamwise advection by the mean profile, diffusion under molecular viscosity and the off-diagonal production term:
\begin{equation} \boldsymbol{A}_{\hat{\boldsymbol{u}}} ={-}\boldsymbol{I} ik_x U + \boldsymbol{A}_{\hat{\boldsymbol{u}},\nu} + \boldsymbol{A}_{\hat{\boldsymbol{u}},P}. \end{equation}Rescaling into the self-similar variables, this becomes
\begin{equation} \boldsymbol{A}_{\hat{\boldsymbol{u}}}(k_x,y,k_z) ={-}\boldsymbol{I} ik_x U + k_z^2\tilde{\boldsymbol{A}}_{\hat{\boldsymbol{u}},\nu}(\xi_x,\xi_z) + k_z\tilde{\boldsymbol{A}}_{\hat{\boldsymbol{u}},P}(\xi_x,\xi_z), \end{equation}
where 
$(\tilde {\cdot })$ is an identical operator, with the differential operators replaced by self-similar counterparts (see Appendix C). Similarly, the gradient acting on the pressure and the eddy viscosity diffusion model can be rescaled as
\begin{equation} \boldsymbol{G} = k_z\tilde{\boldsymbol{G}}(\xi_x,\xi_z) \quad \textrm{and} \quad \mathcal{L}_{\nu_t} = k_z\tilde{\mathcal{L}}_{\nu_t}(\xi_x,\xi_z), \end{equation}
where 
$\nu _t \sim y$ is used for the eddy viscosity diffusion model. The fluctuating equation can now be rewritten as
\begin{equation} \frac{\partial \hat{\boldsymbol{u}}}{\partial t} ={-}\boldsymbol{I} ik_x U + k_z[(\tilde{\boldsymbol{A}}_{\hat{\boldsymbol{u}},P} + \tilde{\mathcal{L}}_{\nu_t} )\hat{\boldsymbol{u}} + \tilde{\boldsymbol{G}}\hat{p}] + k_z^2\tilde{\boldsymbol{A}}_{\hat{\boldsymbol{u}},\nu}\hat{\boldsymbol{u}}+\hat{\boldsymbol{f}}. \end{equation}
In the logarithmic layer, the effects of viscous diffusion can be neglected as 
$Re_\tau \rightarrow \infty$. With rescaled states and forcing below
\begin{equation} \tilde{\boldsymbol{u}}=k_z\hat{\boldsymbol{u}}, \quad \tilde{p}=k_z\hat{p}, \quad \tilde{\boldsymbol{f}}=\hat{\boldsymbol{f}}, \end{equation}the evolution equation for the rescaled state is finally given by
\begin{equation} \frac{\partial \tilde{\boldsymbol{u}}}{\partial \tau} ={-}\boldsymbol{I} ik_x U +(\tilde{\boldsymbol{A}}_{\hat{\boldsymbol{u}},P} + \tilde{\mathcal{L}}_{\nu_t} )\tilde{\boldsymbol{u}} + \tilde{G}\tilde{p} + \tilde{\boldsymbol{f}}, \end{equation}
where time is rescaled with 
$\mathcal {\tau } = k_zt$. Note that 
$\boldsymbol{I} ik_x U$ here does not remain self-similar due to the mean velocity, 
$U(y) \sim \ln y$. However, since this is anti-self-adjoint, it does not enter the spectral budget.
Following the procedure introduced in § 2.2, the Lyapunov-like equation, which determines the energy budget, is now given in terms of the self-similar variables as
\begin{equation} \lim_{t\to\infty}\langle \tilde{\mathcal{L}}_{b}\hat{\boldsymbol{m}},\hat{\boldsymbol{n}} \rangle = 0, \end{equation}where
\begin{align} \tilde{\mathcal{L}}_{b}& \equiv (\tilde{\boldsymbol{A}}_{\tilde{\boldsymbol{u}},P}+\tilde{\mathcal{L}}_{\nu_t})R_{\tilde{\boldsymbol{u}}\tilde{\boldsymbol{u}}}^{\infty} + R_{\tilde{\boldsymbol{u}}\tilde{\boldsymbol{u}}}^{\infty}(\tilde{\boldsymbol{A}}_{\tilde{\boldsymbol{u}},P} + \tilde{\mathcal{L}}_{\nu_t})^{{{\dagger}}}-\tilde{\boldsymbol{G}} R_{\tilde{\boldsymbol{u}}\tilde{p}}^{\infty} - R_{\tilde{p}\tilde{\boldsymbol{u}}}^{\infty} \tilde{\boldsymbol{G}}^{{{\dagger}}}\nonumber\\ &\quad + ( R_{\tilde{\boldsymbol{u}}\tilde{\boldsymbol{f}}}^{\infty} + R_{\tilde{\boldsymbol{f}}\tilde{\boldsymbol{u}}}^{\infty}) \end{align}and
\begin{equation} R_{\tilde{\boldsymbol{u}}\tilde{\boldsymbol{u}}}^{\infty}=k_z^2 R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty}, \quad R_{\tilde{\boldsymbol{u}}\tilde{p}}^{\infty}= k_z^2 R_{\hat{\boldsymbol{u}}\hat{p}}^{\infty}, \quad R_{\tilde{\boldsymbol{u}}\tilde{\boldsymbol{f}}}^{\infty}= k_z R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{f}}}^{\infty}. \end{equation}Furthermore, following Hwang & Lee (Reference Hwang and Lee2020), taking the trace and considering the statistics or spectra, if the pressure transport is assumed small at the considered integral length scales, the balance of the statistics becomes
\begin{equation} {-}2k_z(\mathcal{D}U)g(\xi_x,\xi_z) + \tilde{{\mathsf{E}}}^{\nu_t}_{ii} + \tilde{{\mathsf{E}}}_{ii}^{f} \approx 0, \end{equation}which rearranges to
\begin{equation} \mathcal{D}U = \frac{1}{y}h(\xi_x,\xi_z), \end{equation}where
\begin{equation} h(\xi_x,\xi_z)=\frac{\tilde{{\mathsf{E}}}^{\nu_t}_{ii} + \tilde{{\mathsf{E}}}_{ii}^{f}}{2\xi_z}. \end{equation}Equation (5.17a) can only hold if both sides of it are a constant, resulting in a logarithmic mean profile. Hence, the eddy viscosity dissipation, as well as the self-similar forcing transport term, is essential in making the linear model consistent with a logarithmic mean profile, as is done in the DQLA (Holford et al. Reference Holford, Lee and Hwang2024).
5.3. Colour of turbulence and pressure
It is finally worth pointing out that concerning the modelling of pressure, considering a coloured-in-time modification to the linear operator for the nonlinear term is the only way to have a meaningful model for the slow pressure. Note that any part of the nonlinearity replaced by an irrotational forcing goes straight into the pressure and cannot affect the velocity: the pressure can simply be redefined relative to this term. This is most clearly seen through the pressure Poisson equation:
\begin{equation} \nabla^2 p ={-} 2\frac{\textrm{d}U}{\textrm{d}\kern0.05em y}\frac{\partial v'}{\partial x} + \boldsymbol{\nabla} \boldsymbol{\cdot} \mathcal{N}. \end{equation}
Consider a coloured nonlinear term 
$\mathcal{N} = \mathcal{L}_c\boldsymbol{u}' + \boldsymbol{f}$, where 
$\mathcal {L}_c$ is the operator that adds the colour and 
$\boldsymbol{f}$ is a white-in-time forcing (Zare et al. Reference Zare, Jovanović and Georgiou2017). Since the slow pressure is defined as the component of the pressure driven by 
$\boldsymbol{\nabla} \boldsymbol{\cdot} \mathcal {N}$, the operator 
$\mathcal {L}_c$ can now be used to model the intercomponent transfer of energy mediated by it: i.e. the slow pressure–strain effects. Setting the irrotational forcing to zero, the rapid and slow pressure terms can be determined by
\begin{equation} (\boldsymbol{I} - \boldsymbol{P})\boldsymbol{A}_{\hat{\boldsymbol{u}}} R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty} + R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty}\boldsymbol{A}_{\hat{\boldsymbol{u}}}^{{{\dagger}}}(\boldsymbol{I} - \boldsymbol{P}) = \boldsymbol{G R}_{\hat{\boldsymbol{u}}\hat{p}_r}^\infty + R_{\hat{p}_r\hat{\boldsymbol{u}}}^\infty \boldsymbol{G}^{{{\dagger}}} \end{equation}and
\begin{equation} (\boldsymbol{I} - \boldsymbol{P})\mathcal{L}_c R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty} + R_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty}\mathcal{L}_c^{{{\dagger}}}(\boldsymbol{I} - \boldsymbol{P}) = \boldsymbol{G R}_{\hat{\boldsymbol{u}}\hat{p}_s}^\infty + R_{\hat{p}_s\hat{\boldsymbol{u}}}^\infty \boldsymbol{G}^{{{\dagger}}}, \end{equation}
where 
$\hat {p}_r$ and 
$\hat {p}_s$ are the rapid and slow pressure, respectively (see Appendix B for the derivation). Hence, without the colouring operator, 
$\mathcal {L}_c$, the modelling of slow pressure is not possible.
As it currently stands, an isotropic eddy viscosity diffusion model does not possess the physics to even phenomenologically mimic this. Although not shown, it is found that the pressure–strain arising from the eddy viscosity diffusion model is negligible when compared with the rapid pressure–strain, with the spectra obtained from (5.19a). In other words, the rapid pressure–strain spectra of the DQLA model are virtually indistinguishable from figure 9. Although this may appear to suggest that the pressure–strain arising in the DQLA is not compatible with that from the DNS, decomposing the pressure–strain from the DNS into its rapid and slow components reveals that the pressure–strain from the DQLA does, in fact, resemble that from the DNS.
Figure 16 shows the spanwise wavenumber spectra of the rapid pressure–strain obtained from DNS. It is seen that the signs of rapid pressure–strain spectra for each component in the DNS are consistent with those of the DQLA (compare figures 16 and 9a–c). Most noticeably, the wall-normal component changes sign, indicating that the role of the rapid pressure in the DNS is in line with that in the DQLA. That is, both the streamwise and wall-normal components mainly transfer energy from the source terms to the spanwise components. Moreover, the near-wall ‘splat’ region, or ‘antisplat’ in this case, of the DNS (figure 16b) is also consistent with the DQLA pressure–strain (compare the near-wall red region in figure 16b with that in figure 9b) – the wall-normal component transfers energy to the streamwise component, as opposed to the opposite sign of transfer from the wall-normal to the spanwise component in the slow-pressure case (see figure 9e). The magnitude of the streamwise pressure–strain component is also more consistent with the DQLA, with a majority of the production being redistributed across the components through the slow pressure–strain. The DQLA model with the use of eddy viscosity can only model this rapid pressure–strain. A more physics-rich colouring of the turbulence model for the nonlinear term is required to model slow pressure–strain, a feature responsible for a bulk of the redistribution of production across components.
The premultiplied spanwise wavenumber spectra for the rapid pressure–strain from DNS for the (a) streamwise, (b) wall-normal and (c) spanwise components at 
$Re_\tau \approx 5200$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in 10 % increments. The solid lines are 
$y = 0.01\lambda _x$ and the dashed lines are 
$y = 0.35\lambda _x$.
6. Concluding remarks
In this study, the spectral budget of the DQLA in Holford et al. (Reference Holford, Lee and Hwang2024) was computed and compared with that of DNS. Turbulent fluctuations in the DQLA model are described by the linearised Navier–Stokes equations, with an eddy viscosity diffusion model and a stochastic forcing term to replicate the nonlinear term from the full Navier–Stokes equations. The self-consistent nature of the DQLA ensures that the shear stress generated by the linearised model is consistent with the Reynolds shear stress present in the mean momentum equations. This allows for a fair comparison of the linear fluctuation models’ capabilities across an entire range of length scales. The main findings of this comparison can be summarised as follows:
(a) The budget of the DQLA compares most favourably for the streamwise-elongated structures with strong agreements in production and viscous transport. This is due to the linear model in the DQLA that contains all the necessary physics to replicate such features in combination with the self-consistent nature of the DQLA.
(b) The dissipation under molecular viscosity in the DQLA is largely negligible above the near-wall region. Instead, energy is dissipated locally through the eddy viscosity diffusion term. Overall, the linear model and DQLA framework capture the streamwise-elongated features of turbulence but do not return to isotropic dissipation across the components at smaller length scales, lacking any Kolmogorov scaling features in dissipation, as expected.
(c) By computing the pressure–strain, the key physical process missing from the linear model used in the DQLA becomes apparent: the streak instability and the related nonlinear mechanisms for the generation of streamwise vortical structures. These processes mainly facilitate an efficient transfer of energy from the streamwise component to the other components. This is not available in the linear fluctuation model for the DQLA. In the DQLA, the forcing term, dominated by the wall-normal component, and turbulence production in the streamwise component form the source terms for turbulent kinetic energy. The pressure–strain mainly redistributes this energy into the spanwise component, with some transferred to the streamwise component, resulting in the region of negative production.
(d) There was good agreement between the wall-normal transport by the eddy viscosity model in the DQLA and that of the wall-normal turbulent transport in DNS. The use of the derivative of the eddy viscosity term to replicate the wall-attached features in the wall-parallel components is reaffirmed following Symon et al. (Reference Symon, Madhusudanan, Illingworth and Marusic2023).
(e) Dissipation through the eddy viscosity model provides an essential role in energy removal at integral scales. While the eddy viscosity diffusion model is not energy conservative, its dissipation also keeps the local energy balance self-similar at a considered wavenumber pair. This enables the linear model used in the DQLA to be theoretically consistent with the key properties of the mean-momentum equation in the logarithmic layer, allowing for the self-consistent coupling with the mean equation, as numerically demonstrated by previous studies (Hwang & Eckhardt Reference Hwang and Eckhardt2020; Skouloudis & Hwang Reference Skouloudis and Hwang2021; Holford et al. Reference Holford, Lee and Hwang2024).
Acknowledgements
J.J.H. is supported by Doctoral Training Partnership (DTP) scholarship from the Engineering and Physical Science Research Council (EPSRC) in the UK. Y.H. gratefully acknowledges the support of the Leverhulme Trust (RPG-123-2019) and EPSRC (EP/T009365/1).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Details of the budget formulation
The probabilistic setting for which the Lyapunov-like equations are derived follows the adjoint formulation in Balakrishnan (Reference Balakrishnan1981) (see also Hoepffner et al. Reference Hoepffner, Chevalier, Bewley and Henningson2005; Zare et al. Reference Zare, Jovanović and Georgiou2017; Jovanović Reference Jovanović2021; Holford et al. Reference Holford, Lee and Hwang2024). The covariance of two zero-mean random variables, 
$\xi _1$, 
$\xi _2$, is given by
\begin{equation} \textrm{Cov}_{ \xi_1, \xi_2}(\phi_1, \phi_2) = \mathbb{E}[\langle \xi_1, \phi_1 \rangle_1\langle \xi_2, \phi_2\rangle _2^*], \end{equation}
where 
$(\cdot )^*$ denotes the complex conjugate, 
$\mathbb {E}[\cdot ]$ denotes the expectation operator, 
$\langle \cdot, \cdot \rangle _1$ and 
$\langle \cdot, \cdot \rangle _2$ are the inner products for the corresponding random variables and 
$\phi _1$ and 
$\phi _2$ are arbitrary elements from the same vector space. Given that this covariance function is bilinear, it can be determined through a covariance operator 
$R_{ \xi _1 \xi _2}$, satisfying
\begin{equation} \textrm{Cov}_{ \xi_1, \xi_2}(\phi_1, \phi_2) = \langle R_{ \xi_1 \xi_2} \phi_1, \phi_2 \rangle_2, \end{equation}
where the first of the subscripts denotes the domain and the latter the range of the linear operator. For a budget equation, where most terms are linear transformations of the velocity state, the covariance operator for the various terms is similarly determined through linear transformations and their adjoints. For linear transforms of random variables, say, 
$\tilde {\xi }_{{1}} = \mathcal {L}_{1\tilde {1}}\xi _1$, 
$\tilde {\xi }_{{2}} = \mathcal {L}_{2\tilde {2}}\xi _2$, the covariance operator becomes
\begin{equation} R_{ \tilde{\xi}_{{1}} \tilde{\xi}_{{2}}} = \mathcal{L}_{2\tilde{2}} R_{ \xi_1 \xi_2} \mathcal{L}_{1\tilde{1}}^{{{\dagger}}}, \end{equation}
where 
$(\cdot )^{{\dagger} }$ denotes the adjoint of the linear operator, defined as follows:
\begin{equation} \langle \mathcal{L}_{1\tilde{1}}^{{{\dagger}}}\tilde{\phi}_{{1}}, \phi_{1} \rangle_1 = \langle \mathcal{L}_{1\tilde{1}}\phi_{{1}}, \tilde{\phi}_{{1}} \rangle_{\tilde{1}}, \end{equation}
where 
$\tilde {\phi }_{{1}}$ are elements in the same vector space as 
$\tilde {\xi }_{{1}}$.
To transform the budget equation for 
$\hat {\boldsymbol{u}}$ to that for 
$\hat{\boldsymbol{q}}$, the linear relationship between the states can be used, as well as the operations which remove pressure from the momentum equations. The operations describing the elimination of pressure from the momentum equations are described by the first row, 
$\boldsymbol{B}_{\hat {\boldsymbol{u }}\hat {v}}$, of the following 
$\boldsymbol{B}$ operator:
\begin{equation} \boldsymbol{B} = \begin{bmatrix} \boldsymbol{B}_{\hat{\boldsymbol{u }}\hat{v}} \\ \boldsymbol{B}_{\hat{\boldsymbol{u }}\hat{\eta}} \end{bmatrix} = \begin{bmatrix} -\varDelta^{{-}1} & 0 \\ 0 & \boldsymbol{I} \end{bmatrix}\begin{bmatrix} {\rm i}k_x\mathcal{D} & k^2 & {\rm i}k_z\mathcal{D} \\ {\rm i}k_z & 0 & -{\rm i}k_x \end{bmatrix}. \end{equation}
Here, the divergence is taken, the wall-normal derivative is then applied and finally the divergence of the wall-normal momentum equation is subtracted. The inverse Laplacian is then applied to get the evolution of 
$\hat {v}$. The second row, 
$\boldsymbol{B}_{\hat {\boldsymbol{u }}\hat {\eta }}$, transforms 
$\hat {\boldsymbol{u}}$ into the wall-normal vorticity state 
$\hat {\eta }$, by definition, with the vorticity evolution not involving pressure. To recover the 
$\hat {\boldsymbol{u}}$ state from the 
$\hat{\boldsymbol{q}}$ state, continuity is enforced, resulting in the following:
\begin{equation} \hat{\boldsymbol{u}} = \boldsymbol{C} \hat{\boldsymbol{q}}, \end{equation}where
\begin{equation} \boldsymbol{C} = \begin{bmatrix} C_{\hat{\boldsymbol{q}}\hat{u}} \\ C_{\hat{\boldsymbol{q}}\hat{v}} \\ C_{\hat{\boldsymbol{q}}\hat{w}} \end{bmatrix} = \frac{1}{k^2}\begin{bmatrix} {\rm i}k_x\mathcal{D} & -{\rm i}k_z \\ k^2 & 0 \\ {\rm i}k_z\mathcal{D} & {\rm i}k_x \end{bmatrix}. \end{equation}
This 
$\boldsymbol{C}$ operator is also used to prescribe an inner product for the 
$\hat{\boldsymbol{q}}$ state so that the inner product induces a norm proportional to the turbulent kinetic energy with
\begin{equation} \|\hat{\boldsymbol{u}} \|^2 = \langle \boldsymbol{C} \hat{\boldsymbol{q}}, \boldsymbol{C} \hat{\boldsymbol{q}} \rangle = \langle \boldsymbol{M} \hat{\boldsymbol{q}}, \hat{\boldsymbol{q}} \rangle = \langle \hat{\boldsymbol{q}}, \hat{\boldsymbol{q}} \rangle_e = \|\hat{\boldsymbol{q}}\|^2_e, \end{equation}
where 
$\langle \cdot, \cdot \rangle$ is the unweighted standard vector inner product of the corresponding vector space and
\begin{equation} \boldsymbol{M} = \frac{1}{k^2}\begin{bmatrix} -\varDelta & 0 \\ 0 & \boldsymbol{I} \end{bmatrix} \end{equation}
is a positive definite matrix operator to weight the 
$\hat{\boldsymbol{q}}$ vector space. With this inner product, using the adjoint definition (A3), 
$\boldsymbol{B}^{{\dagger} } = \boldsymbol{C}$ (Jovanović & Bamieh Reference Jovanović and Bamieh2005). The equation for the steady-state covariance operator of 
$\hat{\boldsymbol{q}}$ can then be determined from (2.8) by applying 
$\boldsymbol{B}$ and using (A2) and (A5a) with
\begin{equation} \boldsymbol{B} \mathcal{L}_b^{\infty}\boldsymbol{B}^{{{\dagger}}} = \boldsymbol{A}_{\hat{\boldsymbol{q}}} R_{\hat{\boldsymbol{q}}\hat{\boldsymbol{q}}}^{\infty} + R_{\hat{\boldsymbol{q}}\hat{\boldsymbol{q}}}^{\infty} \boldsymbol{A}_{\hat{\boldsymbol{q}}}^{{{\dagger}}} + \boldsymbol{BR}_{\hat{\boldsymbol{q}}\hat{\boldsymbol{\mathcal{N}}}}^\infty + R_{\hat{\boldsymbol{\mathcal{N}}}\hat{\boldsymbol{q}}}^\infty \boldsymbol{B}^{{{\dagger}}} = 0. \end{equation}Here,
\begin{equation} \boldsymbol{B} \boldsymbol{A}_{\hat{\boldsymbol{u}}}\boldsymbol{C} \hat{\boldsymbol{q}} = \boldsymbol{A}_{\hat{\boldsymbol{q}}}\hat{\boldsymbol{q}} = \begin{bmatrix} \varDelta^{{-}1}\mathcal{L}_{OS} & 0 \\ -{\rm i}k_z(\mathcal{D}U) & \mathcal{L}_{SQ} \end{bmatrix} \begin{pmatrix} \hat{v} \\ \hat{\eta} \end{pmatrix} \end{equation}is the Orr–Sommerfeld–Squire system, where
\begin{equation} \mathcal{L}_{OS} ={-}{\rm i}k_x(U\varDelta - (\mathcal{D}^2U))+\frac{1}{Re}\varDelta^2 \end{equation}
is the Orr–Sommerfeld operator, where 
$\varDelta ^2 = \varDelta \varDelta$ with the additional boundary conditions 
$\mathcal {D}\hat {v} = 0$ at 
$y=0$ and 
$y=2h$ to satisfy continuity at the wall.
Appendix B. Computation of some cross-correlations
Since the pressure does not explicitly appear in the present formulation based on the Orr–Sommerfeld–Squire system, here the computation method for the forcing–velocity and the pressure–velocity correlations is described. Note that only the solenoidal part of the forcing can affect the velocity, and the 
$\boldsymbol{B}$ operator that acts on the forcing in the Orr–Sommerfeld–Squire system removes the irrotational part of the forcing. From this observation, combining the 
$\boldsymbol{B}$ and 
$\boldsymbol{C}$ operators together into a single expression defines a projection operator,
\begin{equation} \boldsymbol{P} = \boldsymbol{C} \boldsymbol{B}, \end{equation}
eliminating the irrotational field, projecting a vector onto its solenoidal component. Moreover, 
$\boldsymbol{C}^{{\dagger} } = \boldsymbol{B}$ (Jovanović & Bamieh Reference Jovanović and Bamieh2005), where the inner product on the wall-normal velocity vorticity state is given by (A6a). Therefore, 
$\boldsymbol{P}$ is an orthogonal projection operator, ensuring a unique orthogonal decomposition of a vector into its solenoidal and irrotational parts, provided the operator has suitably defined boundary conditions, i.e. no-slip for the velocity state. Therefore, for a three-dimensional vector 
$\hat {\boldsymbol{m}}$, which is sufficiently smooth to allow the following decomposition, 
$\hat {\boldsymbol{m}} = \hat {\boldsymbol{m}}^s + \hat {\boldsymbol{m}}^i$, where 
$\boldsymbol{\nabla} \boldsymbol {\cdot } \boldsymbol m^s = 0$ and 
$\boldsymbol \nabla \times \boldsymbol m^i = 0$, then 
$\boldsymbol m^s$ and 
$\boldsymbol m^i$ can be determined from
$$\begin{gather} \hat{\boldsymbol{m}}^s = \boldsymbol P\hat{\boldsymbol{m}}, \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{m}}^i = (\boldsymbol I - \boldsymbol P)\hat{\boldsymbol{m}}, \end{gather}$$
where 
$\boldsymbol I$ is the identity operator of the corresponding space. Note that this decomposition is referred to as the Helmholtz decomposition.
Application of the Helmholtz decomposition to the linearised fluctuation equations for the DQLA subsequently yields
$$\begin{gather} \frac{\partial{\hat{\boldsymbol{u}}}}{\partial t}=\boldsymbol P (A_{\hat{\boldsymbol{u}}} + \mathcal{L}_{\nu_t})\hat{\boldsymbol{u}}+\hat{\boldsymbol f}^s, \end{gather}$$
$$\begin{gather}\boldsymbol 0 = (\boldsymbol I - \boldsymbol P)(A_{\hat{\boldsymbol{u}}} + \mathcal{L}_{\nu_t})\hat{\boldsymbol{u}} -\boldsymbol G \hat{p}+\hat{\boldsymbol f}^i, \end{gather}$$
where 
$\mathcal {L}_{\nu _t}$ is the eddy-viscosity-related part of 
$\mathcal {N}$ in (3.1). Note that the pressure satisfies the following Poisson equation:
\begin{equation} \hat{\boldsymbol \nabla} \boldsymbol{\cdot} (\boldsymbol G \hat{p}) ={-}2{\rm i}k_x\hat{v}\frac{\textrm{d} U}{\textrm{d}\kern0.05em y} + \hat{\boldsymbol \nabla} \boldsymbol{\cdot} ((\boldsymbol I - \boldsymbol P)\mathcal{L}_{\nu_t} \hat{\boldsymbol{u}}) + \hat{\boldsymbol \nabla} \boldsymbol{\cdot} \hat{\boldsymbol f}^i, \end{equation}
with 
$\hat {\boldsymbol \nabla }=(i k_x, \mathrm {d}/\mathrm {d}\kern 0.05em y, {\rm i} k_z)$. It is useful to further decompose the pressure such that 
$\hat {p}=\hat {p}^u+\hat {p}^f$, where 
$\hat {p}^u$ and 
$\hat {p}^f$ are from the first two terms and the last term of the right-hand side of (B3c), respectively. Since 
$\hat {\boldsymbol f}^i=\hat {\boldsymbol \nabla } \hat {\phi }$ by definition, where 
$\hat {\phi }$ is a scalar potential for 
$\hat {\boldsymbol f}^i$, it becomes evident that 
$\hat {p}^f=\hat {\phi }$. This suggests that 
$\hat {p}^f$ is only directly correlated with the irrotational part of the forcing, which does not influence the velocity field in the DQLA. Therefore, 
$\hat {p}^f$ and 
$\hat {\boldsymbol{u}}$ must not be correlated with each other, indicating that the correlations between the pressure and irrotational part of the forcing in the DQLA must vanish, i.e.
\begin{equation} R_{\hat{p}\hat{\boldsymbol{u}}}=R_{\hat{p}^u\hat{\boldsymbol{u }}} \ \textrm{or, equivalently,} \varPhi_{\hat{p}\hat{\boldsymbol{u}}}=\varPhi_{\hat{p}^u\hat{\boldsymbol{u }}}, \end{equation}
where 
$\varPhi _{\hat {p}\hat {\boldsymbol{u }}}=\mathbb {E}[\kern0.7pt \hat{p}(y)\hat {\boldsymbol{u }}^H(y')]$. Similarly, the correlations between 
$\hat {\boldsymbol{u}}$ and 
$\hat {f}_i$ must vanish, indicating
\begin{equation} R_{\hat{\boldsymbol f}\hat{\boldsymbol{u }}}=R_{\hat{\boldsymbol f}^s\hat{\boldsymbol{u}}}\ \textrm{or, equivalently,} \varPhi_{\hat{\boldsymbol f}\hat{\boldsymbol{u }}}=\varPhi_{\hat{\boldsymbol f}^s\hat{\boldsymbol{u}}}, \end{equation}
where 
$\varPhi _{\hat {\boldsymbol f }\hat {\boldsymbol{u}}}=\mathbb {E}[\,\hat {\boldsymbol f }(y)\hat {\boldsymbol{u}}^H(y')]$.
Now, given 
$\varPhi ^{\infty }_{\hat {\boldsymbol{u}}\hat {\boldsymbol{u}}}$ available from the DQLA model, the forcing–velocity and pressure–velocity correlations can be determined. First, combining with (B4a), right multiplication of (B3a) by 
$\boldsymbol u^H$ and subsequent application of the expectation operation give
\begin{equation} \varPhi_{\hat{\boldsymbol f}\hat{\boldsymbol{u }}}={-}\boldsymbol P(A_{\hat{\boldsymbol{u}}} + \mathcal{L}_{\nu_t}) \varPhi_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty}. \end{equation} Similarly, after right multiplication of (B3a) by 
$\boldsymbol u^H$, taking divergence and the expectation operator to (B3b) yields
\begin{equation} \hat{\nabla}^2 \varPhi_{\hat{p}\hat{\boldsymbol{u}}}^\infty=\hat{\boldsymbol \nabla}\boldsymbol{\cdot} (\boldsymbol I - \boldsymbol P)(A_{\hat{\boldsymbol{u}}} + \mathcal{L}_{\nu_t})\varPhi_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}^{\infty}, \end{equation}using (B4b). The pressure–velocity correlation is finally obtained by inverting the Laplacian above subject to Neumann boundary conditions.
Appendix C. Definition of the rescaled self-similar operators
To rescale the linear operators in terms of self-similar variables (5.5), using 
$k_z$ as a reference length scale, the differential operators become
\begin{equation} {\rm i}k_x = {\rm i}k_z\frac{\xi_x}{\xi_z},\quad \mathcal{D} = k_z \frac{\textrm{d}}{\textrm{d}\xi_z} = k_z\tilde{D}, \end{equation}
where 
$\tilde {D} = {\mathrm {d}}/{\mathrm {d}\xi _z}$ and 
${\rm i}k_z$ is a constant. The gradient can then be expressed as
\begin{equation} \boldsymbol G(k_x,y,k_z) = k_z \tilde{\boldsymbol G}(\xi_x,\xi_z), \end{equation}where
\begin{equation} \tilde{\boldsymbol G}(\xi_x,\xi_z) = \begin{bmatrix} {\rm i}\dfrac{\xi_x}{\xi_z} & \tilde{\mathcal{D}} & {\rm i} \end{bmatrix}^{\rm T}. \end{equation}The Laplacian is given by
\begin{equation} \varDelta(k_x,y,k_z) = k_z^2\tilde{\varDelta}(\xi_x,\xi_z), \end{equation}where
\begin{equation} \tilde{\varDelta}(\xi_x,\xi_z) = \tilde{\mathcal{D}}^2 - \left(\frac{\xi_x^2}{\xi_z^2}+1\right). \end{equation}Similarly, the mean shear and the eddy viscosity can be rescaled approximately in the logarithmic layer. Using the logarithmic mean profile, the mean shear is rescaled as
\begin{equation} (\mathcal{D}U)(y) \approx \frac{1}{\kappa y} = k_z\widetilde{\mathcal{D}U}(\xi_z), \end{equation}
with 
$\widetilde {\mathcal {D}U}(\xi _z)=1/(\kappa \xi _z)$. The rescaling of the eddy viscosity relies upon the mean momentum equation (5.3) in the logarithmic layer. Using 
$\nu _t(\mathcal {D}U)(y) = u_\tau ^2$,
\begin{equation} \nu_t = u_\tau^2\kappa y = k_z^{{-}1}\tilde{\nu}_t \end{equation}is obtained, where
\begin{equation} \tilde{\nu}_t = u_\tau^2 \kappa \xi_z \end{equation}and
\begin{equation} \mathcal{D}\nu_t(y) = u_\tau^2\kappa = \widetilde{\mathcal{D}\nu_t}(\xi_z). \end{equation} Substituting these into the 
$\boldsymbol{A}_{\hat {\boldsymbol{u}}}$ decomposition (5.9) gives
$$\begin{gather} \boldsymbol{A}_{\hat{\boldsymbol{u}},P}(k_x,y,k_z) = k_z\begin{bmatrix} 0 & \widetilde{\mathcal{D}U} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} = k_z \widetilde{\boldsymbol A}_{\hat{\boldsymbol{u}},P}(\xi_z), \end{gather}$$
$$\begin{gather}\boldsymbol{A}_{\hat{\boldsymbol{u}}, \nu}(k_x,y,k_z)= k_z^2\begin{bmatrix} \tilde{\varDelta} & 0 & 0 \\ 0 & \tilde{\varDelta} & 0 \\ 0 & 0 & \tilde{\varDelta} \end{bmatrix} = k_z^2 \widetilde{\boldsymbol A}_{\hat{\boldsymbol{u}},\nu}(\xi_x,\xi_z), \end{gather}$$and the eddy viscosity operator becomes
\begin{equation} \mathcal{L}_{\nu_t}(k_x,y,k_z) = k_z\begin{bmatrix} \widetilde{\mathcal{D}\nu_t}\tilde{\mathcal{D}} + \tilde{\nu}_t\tilde{\varDelta} & {\rm i}\dfrac{\xi_x}{\xi_z}\widetilde{\mathcal{D}\nu_t} & 0\\ 0 & 2\widetilde{\mathcal{D}\nu_t}\tilde{\mathcal{D}} + \tilde{\nu}_t\tilde{\varDelta} & 0 \\ 0 & {\rm i} \widetilde{\mathcal{D}\nu_t} & \widetilde{\mathcal{D}\nu_t}\tilde{\mathcal{D}} + \tilde{\nu}_t\tilde{\varDelta} \end{bmatrix} = k_z \tilde{\mathcal{L}}_{\nu_t}(\xi_x,\xi_z). \end{equation}Appendix D. Effects of $N_{POD}$ on turbulent transport and dissipation in DQLA
To further justify the use of 
$N_{POD} = 2$, figure 17 shows the spanwise wavenumber spectra of total turbulent transport and dissipation from the DQLA. For 
$N_{POD} = 1530$, it can be seen that the main role of the nonlinear model term in the DQLA (figure 17a) is providing energy at small length scales (
$\lambda _{z/h}\approx 0.5$), as the turbulent transport spectra are mostly positive especially around the small spanwise length scales. In the dissipation spectra (figure 17b), it can be seen that most of this energy is dissipated locally. However, this feature appears across the entire channel and is particularly strong at the channel centre. The strong positive turbulent transport in the channel centre and the corresponding dissipation are non-physical, and this is evident when comparing with the interscale transport spectra of DNS (figure 15d–f). Upon decreasing 
$N_{POD}$, this non-physical artefact remains for 
$N_{POD} = 8$ (figure 17c), with the dissipation spectra also extending to the channel centre (figure 17d). For the 
$N_{POD} = 2$ case, the positive turbulent transport near the channel centre (3.1) is almost absent, having little effect on the dissipation (figure 17e, f).
The premultiplied spanwise wavenumber spectra for (a,c,e) total turbulent transport and for (b,d,f) dissipation: (a,b) 
$N_{POD} =1530$; (c,d) 
$N_{POD} = 8$; (e, f) 
$N_{POD} = 2$. Contour levels are separated by 1 % of the maximum value up to 10 % of the maximum, and then in 10 % increments.
Open access
