2.1. Justifying the single-scale analysis

The present work develops a new model for predicting phase profiles of scale interactions based on the analytical foundations of Jacobi et al. (Reference Jacobi, Chung, Duvvuri and McKeon2021). In order to develop this new model, it is necessary to first revisit those analytical foundations. In particular, the traditional approach to quantifying scale interactions involves filtering large- and small-scale velocity signals and computing correlations between them, whereas the analytical resolvent approach considers correlations of velocity and Reynolds stress fluctuations at a single, isolated scale, without filtering to separate between large- and small-scale motions. At first glance, these two techniques seem quite different. The single-scale approach assumed that isolated, VLSM-scale velocity and stress fluctuations are the dominant contributors to inter-scale interactions and thus capture the relevant trends obtained from the traditional filtered correlations.

To justify the assumption of the single-scale approach, we performed a spectral comparison of the correlation coefficient of Mathis et al. (Reference Mathis, Hutchins and Marusic2009), calculated by filtering, and an analogous correlation coefficient based on single-scale velocity and Reynolds stress calculations, without filtering. In this way, we can identify the conditions under which the velocity and stress signals (without filtering) can replicate trends obtained from filtering the velocity signal, and also when the single scale of the velocity and stress signals can capture the integrated effect of a wide range of scales. These conditions will then provide a physical context to the resolvent-based model developed in the following sections.

Throughout the paper, we non-dimensionalize the incompressible Navier–Stokes equations with respect to the friction velocity $u_\tau$ and outer length scale, $h$, corresponding to a half-channel height that is relevant for later comparisons with channel DNS calculations. The Reynolds number is defined as ${Re} \equiv {u_\tau h}/{\nu }$, where $\nu$ is the kinematic viscosity; $U_0$ denotes the centreline velocity and $(x,y,z)$ and $(u,v,w)$ are the streamwise, wall-normal and spanwise coordinates and velocities, respectively.

Mathis et al. (Reference Mathis, Hutchins and Marusic2009) defined their amplitude modulation coefficient in terms of a large-scale signal, $u_L$, obtained by low-pass filtering the instantaneous velocity, $u$, and a small-scale remainder signal, $u_S$, obtained by subtraction, $u_S = u - u_L$. The envelope of the small scales was then calculated by the Hilbert transform in order to describe the large-scale manifestation of the small-scale fluctuations; in our analysis, we calculate a roughly equivalent envelope by squaring the small-scale signal, $u_S^2$, following the approach of Chung & McKeon (Reference Chung and McKeon2010). The resulting large-scale and envelope signals were then Fourier transformed into wavenumber space, denoted by the $\widehat {( {\cdot } )}$, to yield the amplitude modulation coefficient via the cross-correlation theorem

(2.1) \begin{equation} R(y) = \frac{ \mathbb{R} [ \int_{-\infty}^\infty \hat{u}_L^* \widehat{u^2_S} \, {\rm d}k_x ]}{[ \int_{-\infty}^\infty |\hat{u}_L|^2 \, {\rm d}k_x ]^{1/2} [ \int_{-\infty}^\infty |\widehat{u^2_S}|^2 \, {\rm d}k_x]^{1/2}}, \end{equation}

where $^*$ denotes the complex conjugate. Rewriting the complex signals, at each wavenumber $k_x$, in terms of their magnitudes and the phase difference between them, $\Delta \phi (k)$, according to

(2.2)\begin{equation} \hat{u}_L^* \widehat{u^2_S} = |\hat{u}_L| |\widehat{u^2_S}| {\rm e}^{{\rm i} \Delta \phi} = |\hat{u}_L| |\widehat{u^2_S}| \left[\cos{(\Delta \phi)} + {\rm i} \sin{(\Delta \phi)} \right], \end{equation}

we can then rewrite the modulation coefficient as

(2.3)\begin{equation} R(y) = \frac{ \int_{-\infty}^\infty |\hat{u}_L| |\widehat{u^2_S}| \cos{(\Delta \phi)} \, {\rm d}k_x }{[ \int_{-\infty}^\infty |\hat{u}_L|^2 \, {\rm d}k_x ]^{1/2} [ \int_{-\infty}^\infty |\widehat{u^2_S}|^2 \, {\rm d}k_x ]^{1/2}} = \langle \cos{(\Delta \phi)} \rangle_w. \end{equation}

From this integral, the amplitude modulation coefficient can be naturally interpreted as an average phase cosine, $\langle \cos {(\Delta \phi )} \rangle _w$, weighted by the magnitude of the cross-spectral density, $\hat {u}_L^* \widehat {u^2_S}$. Therefore, the amplitude modulation coefficient can be interpreted physically as an average phase difference between the large- and small-scale signals. (Note that the definition of the cross-spectrum here differs in the complex conjugate from that defined in Jacobi & McKeon (Reference Jacobi and McKeon2013) in order to make the phase difference consistent with Jacobi et al. (Reference Jacobi, Chung, Duvvuri and McKeon2021).)

Figure 1(a) shows the premultiplied cross-spectral energy density $|\hat {u}_L^* \widehat {u^2_S}|$ in the inner region of a channel flow DNS at ${Re} = 5200$ by Lee & Moser (Reference Lee and Moser2015), ensemble averaged over 15 300 streamwise/wall-normal velocity-field slices. The cross-spectral energy density represents the weighting factor for the phase cosine and the integral of the normalized cross-spectrum yields the amplitude modulation coefficient, shown as the solid blue line in figure 1(b). Despite the different envelope used here, the general shape of the modulation profile is consistent with previous reports of Mathis et al. (Reference Mathis, Hutchins and Marusic2009) and others, specifically showing a zero-crossing location near the wall, which represents the location of a $-{\rm \pi} /2$ average phase difference between the two signals. Except for the signature of the near-wall cycle, most of the co-spectral energy in the inner region is concentrated in the large-scale motions, for $0.4 \lessapprox k_x \lessapprox 2$, consistent with experiments by Jacobi & McKeon (Reference Jacobi and McKeon2013).

Figure 1. Magnitude of the pre-multiplied cross-spectrum of (a) $| \langle k_x \hat {u}^*_L \widehat {u^2_S} \rangle |$ following Mathis et al. (Reference Mathis, Hutchins and Marusic2009) with a root mean square (r.m.s.) envelope, and (c) $| \langle k_x \hat {u}^* \widehat {u^2} \rangle |$ used in the present study. The two white dashed lines represent the large-scale region $0.4 < k_x < 2$. Integrating the cross-spectrum over all wavenumbers yields the amplitude modulation coefficient, $R(y)$, shown in (b) for both cases: the filtered cross-spectrum of (a) shown in blue, and the unfiltered cross-spectrum of (c) shown in red. The difference between the two coefficient profiles is no greater than the variability due to the choice of filter cutoff explained in Mathis et al. (Reference Mathis, Hutchins and Marusic2009). The coefficient $R(y)$ resulting from integrating over only the large-scale sub-region in (c) is shown as the dashed line, and nearly coincides with the full integral in red, except very near the wall.

Figure 1(c) shows the premultiplied co-spectral energy density for the raw, unfiltered velocity and stress signals, $|\hat {u}^* \widehat {u^2}|$, to compare with the filtered results in (a). Both spectral maps show the same dominant contribution from the VLSMs, although the unfiltered case is more intense. The similarity between the filtered and unfiltered cross-spectral maps at the large scales can be explained by considering their algebraic relationship

(2.4) \begin{equation} \underbrace{\vphantom{\left(\widehat{u^2_S}^* + 2 \widehat{u_L u_S}^*+ \widehat{u^2_L}^* \right)} \hat{u}^* \widehat{u^2}}_\text{unfiltered} = \underbrace{\vphantom{\left(\widehat{u^2_S}^* + 2 \widehat{u_L u_S}^*+ \widehat{u^2_L}^* \right)}\hat{u}^*_L \widehat{u^2_S}}_\text{filtered} + \underbrace{\vphantom{\left(\widehat{u^2_S}^* + 2 \widehat{u_L u_S}^*+ \widehat{u^2_L}^* \right)}2 \hat{u}^*_L \widehat{u_L u_S}}_\text{negligible} + \underbrace{\vphantom{\left(\widehat{u^2_S}^* + 2 \widehat{u_L u_S}^*+ \widehat{u^2_L}^* \right)}\hat{u}^*_L \widehat{u^2_L}}_\text{self-modulation} + \underbrace{\vphantom{\left(\widehat{u^2_S}^* + 2 \widehat{u_L u_S}^*+ \widehat{u^2_L}^* \right)}\hat{u}^*_S(\widehat{u^2_S} + 2 \widehat{u_L u_S}+ \widehat{u^2_L} )}_\text{small-scale contribution}. \end{equation}

The terms that multiply $\hat {u}_S$ contribute energy only at the small scales and can be neglected when considering interactions at small wavenumbers. The remaining terms that explain the discrepancy between the filtered and unfiltered signals are (i) a cross-interaction term $\hat {u}_L^* \widehat {u_L u_S}$ which contributes very little due to the lack of spectral overlap between the large- and small-scale signals, and (ii) a self-modulation term, $\hat {u}^*_L \widehat {u^2_L}$, which represents interactions between large scales and the envelope of other large scales on the same side of the filter cutoff.

The relative energetic balance between the filtered scale-interaction term $\hat {u}^*_L \widehat {u^2_S}$ and the self-interaction term $\hat {u}^*_L \widehat {u^2_L}$ determines the extent to which the unfiltered cross-spectrum mimics the filtered cross-spectrum, and this balance is a result of the choice of filter cutoff. As the filter cutoff increases, the cross-spectral energy at the large scales decreases and Mathis et al. (Reference Mathis, Hutchins and Marusic2009) noted that the magnitude of the modulation effect also decreases. This decrease in cross-spectral energy is due to the fact that the modulation effect primarily represents interactions between VLSMs and scales that are (relatively speaking) smaller, but not strictly small. If the filter-cutoff wavenumber is too high, then the only scales in the envelope signal are too small to interact with the VLSMs, and all of the cross-spectral energy shifts to the self-modulation term. As the filter-cutoff wavenumber is decreased, the energy in the self-modulation also decreases and the filtered and unfiltered scale-interaction terms become more aligned, as shown in detail in Appendix A. As long as the filter-cutoff wavenumber is sufficiently low, the unfiltered cross-spectrum will capture the same interactions between VLSMs and relatively smaller scales as the filtered cross-spectrum.

Therefore, instead of filtering the instantaneous velocity signals to construct a cross-spectrum of $\hat {u}_L^* \widehat {u^2_S}$, we can actually skip the filtering step altogether and construct a cross-spectrum based on the instantaneous velocity and stress signals, $\hat {u}^* \widehat {u^2}$, and the discrepancy between these two representations should be of the same order of magnitude as the variability observed in $\hat {u}_L^* \widehat {u^2_S}$ alone due to changing the filter cutoff (Mathis et al. Reference Mathis, Hutchins and Marusic2009). The resulting amplitude modulation coefficient from the unfiltered signals, shown as the solid red line in figure 1(b), exhibits all of the same qualitative behaviour as the coefficient calculated from the filtered signals, and a quantitatively similar zero-crossing location.

In addition to showing that the velocity/stress cross-spectrum can approximate the filtered cross-spectrum, figure 1 also shows that only a narrow range of wavenumbers contributes significantly to the scale-interaction effect, marked in the region between the two vertical dashed lines. Integrating the cross-spectrum within just this wavenumber range produces an amplitude modulation coefficient in figure 1(b) that is nearly identical to the coefficient resulting from integrating over all wavenumbers. Therefore, we can focus on just isolated VLSM wavenumbers in the resolvent analysis without filtering, and these should represent the dominant contribution to the overall scale-interaction behaviour.

The use of unfiltered velocity and stress signals appears to be consistent with the traditional filter-based analysis of scale interactions, at least for VLSMs interacting with smaller scales; and the assumption that isolated wavenumbers can capture the integrated effect appears to apply except very close to the wall. Figure 1(c) shows that the VLSM contribution to $R(y)$ varies with wall-normal location. Beneath the log layer ($y<0.025$ or $y^+ < 130$), the cross-spectral energy density is broadly distributed across a wide range of wavenumbers, so the single-scale modelling will be less accurate. Special techniques will be developed in § 4 to unwrap the phase difference between the large-scale velocity mode and Reynolds stress fluctuations at VLSM scales in this near-wall region.

2.2. Triple-decomposed equations

Having established that the single-scale representations of velocity and stress can be used as a proxy for studying large- and small-scale turbulent interactions, we now briefly explain the resolvent-based representation for scale interactions used in Jacobi et al. (Reference Jacobi, Chung, Duvvuri and McKeon2021), which was based on earlier work by Reynolds & Hussain (Reference Reynolds and Hussain1972).

A triple decomposition of the instantaneous velocity, $u_i$, and pressure, $p$, fields was then performed, separating the ensemble average values, denoted $\overline {(\;)}$, from a single, isolated coherent scale, denoted $\widetilde {(\;)}$, and the remaining turbulent components, denoted $(\;)'$, according to

(2.5a,b)\begin{equation} u_i = \bar{u}_i + \tilde{u}_i + u_i', \quad p = \bar{p} + \tilde{p} + p'. \end{equation}

This triple decomposition can be accomplished practically by a very narrow, bandpass spectral filter around a single (vector-valued) wavenumber triplet of interest, $\boldsymbol {k_f}=(k_x,k_z,\omega )$, which will reflect the large-scale wavenumber of the VLSMs. The components of $\boldsymbol {k_f}$ are the real wavenumbers in the wall-parallel directions, $k_x$ and $k_z$, and the real angular frequency of the associated coherent structure, $\omega$. Thus $\tilde {u}_i$ represents just a single, isolated velocity mode, and $u_i'$ represents all of the remaining wavenumber triplets in the turbulence signal, $\boldsymbol {k} \neq \boldsymbol {k_f}$.

Substituting the triple decomposition (2.5a,b) into the incompressible Navier–Stokes equations results in a dynamical equation that includes an instantaneous Reynolds stress term, $u_i' u_j'$, which represents content from all wavenumber triplets, including $\boldsymbol {k_f}$, due to the quadratic nonlinearity. Filtering the resulting equation at the isolated scale, $\boldsymbol {k_f}$, results in a dynamical description of just the isolated scale motion, $\tilde {u}_i$

(2.6a,b)\begin{equation} \frac{\partial \tilde{u}_i}{\partial t} + \bar{u}_j\frac{\partial \tilde{u}_i}{\partial x_j} + \tilde{u}_j\frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \tilde{p}}{\partial x_i} - \frac{1}{Re} \frac{\partial^2 \tilde{u}_i}{\partial x_j^2} ={-} \frac{\partial \tilde{r}_{ij}}{\partial x_j} \equiv {\tilde{f}_i}, \quad \frac{\partial \tilde{u}_i}{\partial x_i} = 0, \end{equation}

which depends on the filtered Reynolds stress, $\tilde {r}_{ij} = \widetilde {u_i' u_j'}$.

If we consider just the streamwise component, for illustrative purposes, then the isolated scale, $\tilde {u}$, can be thought of as representing the dominant streamwise VLSM mode. Thus, $\tilde {u}$ can be interpreted as the large-scale, streamwise velocity signal commonly derived in previous studies of inter-scale interactions via low-pass filtering, $u_L$.

Similarly, the instantaneous Reynolds stress, $u'^2$, can be thought of as a variance envelope of the remainder turbulent fluctuations, $u'$, after the VLSM component was removed, as discussed in detail in § 2.1. Filtering the enveloped fluctuations to obtain $\widetilde {u'^2}$ is then equivalent to rectifying the enveloped signal at the same scale as the large-scale signal, $\tilde {u}$. The enveloped, non-VLSM signal, $\widetilde {u'^2}$, differs from the enveloped small-scale signal used in previous studies of scale interactions, $u_S$ or $\mathcal {E}(u_S)$, due to the presence of large-scale (albeit non-VLSM) contributions in $u'$. Nevertheless, $\widetilde {u'^2}$ can still be interpreted as equivalent to the traditional, small-scale envelope assuming that the non-VLSM contributions to the scale interactions are negligible, as shown above in figure 1.

The dynamical equation for the filtered Reynolds stress, $\tilde {r}_{ij}$, was derived in similar fashion

(2.7a)\begin{gather} \frac{\partial \tilde{r}_{ij}}{\partial t} + \bar{u}_k \frac{\partial \tilde{r}_{ij}}{\partial x_k} + \tilde{u}_k \frac{\partial \bar{r}_{ij}}{\partial x_k} + \tilde{r}_{jk} \frac{\partial \bar{u}_i}{\partial x_k} + \tilde{r}_{ik} \frac{\partial \bar{u}_j}{\partial x_k} + \bar{r}_{jk} \frac{\partial \tilde{u}_i}{\partial x_k} + \bar{r}_{ik} \frac{\partial \tilde{u}_j}{\partial x_k} - \frac{1}{Re}\frac{\partial^2 \tilde{r}_{ij}}{\partial x_k^2} = \tilde{g}_{ij}, \end{gather}

(2.7b)\begin{gather} \tilde{g}_{ij} ={-} \frac{\partial }{\partial x_k} {\widetilde{u_i' u_j' u_k'}} - \widetilde{u_j'\frac{\partial {p'}}{\partial x_i}} - \widetilde{u_i'\frac{\partial {p'}}{\partial x_j}} - \frac{2}{Re} \widetilde{\frac{\partial {u}_i'}{\partial x_k}\frac{\partial {u}_j'}{\partial x_k}}, \end{gather}

where the nonlinear terms have been grouped together on the right-hand side, labelled $\tilde {g}_{ij}$. Note that the filtered Reynolds stress dynamical equation depends on the ensemble averaged Reynolds stress profiles, $\bar {r}_{jk} = \overline {u'_j u'_k}$. The channel flow geometry resulted in simplifications of the mean flow

(2.8a,b)\begin{equation} \bar{u} = \bar{u}(y), \quad \{\bar{v}, \bar{w}\} = 0, \end{equation}

and mean Reynolds stresses

(2.9a,b)\begin{equation} \{ \bar{r}_{xx},\bar{r}_{yy},\bar{r}_{zz}, \bar{r}_{xy}\} = \{\bar{r}_{xx},\bar{r}_{yy},\bar{r}_{zz}, \bar{r}_{xy}\}(y), \quad \{\bar{r}_{xz},\bar{r}_{yz}\} = 0, \end{equation}

based on geometrical and statistical symmetries.

2.3. Resolvent formulation

To apply the resolvent analysis to this system of equations, the isolated (large) scale, $\tilde {u}_i$, and filtered Reynolds stress, $\tilde {r}_{ij}$, signals were expressed in the form of normal modes. In addition, the nonlinear forcing of the isolated mode (which appears as the divergence of the filtered Reynolds stress in (2.6)) was also assumed to exhibit normal mode form, according to

(2.10)\begin{equation} {\{\tilde{u}_i, \tilde{r}_{ij}, \tilde{f}_i \}(x,y,z,t) = \{\tilde{U}_i, \tilde{R}_{ij}, \tilde{F}_i\}(y;k_x,k_z,\omega) \exp({\mathrm{i}(k_x x + k_z z - \omega t)}){+ \text{c.c.}}}, \end{equation}

where $\tilde {f}_i$ is the forcing term in (2.6). The complex mode profiles, $\tilde {U}_i$, $\tilde {R}_{ij}$ and $\tilde {F}_i$ can be decomposed in terms of a magnitude and phase as $\tilde {U}_i = |\tilde {U}_i| {\rm e}^{{\rm i} \phi _{U_i}}$, $\tilde {R}_{ij} = |\tilde {R}_{ij}| {\rm e}^{{\rm i} \phi _{R_{ij}}}$ and $\tilde {F}_{i} = |\tilde {F}_{i}| {\rm e}^{{\rm i} \phi _{F_{i}}}$, and the phase difference between large-scale velocity and stress modes is defined as $\Delta \phi = \phi _{\tilde {R}_{xx}}-\phi _{\tilde {U}}$. The $\text {c.c.}$ denotes the complex conjugate of the preceding normal modes.

Substituting the normal mode decomposition (2.10) into the large-scale dynamics (2.6) results in a dynamical equation for each wavenumber $(k_x,k_z,\omega )$. Combining the three components of velocity modes $\tilde {u}_i$ into a matrix form, $\boldsymbol{\mathsf{U}}$, and then re-writing the velocity components in wall-normal velocity/vorticity form, $\boldsymbol{\mathsf{Q}}$ (following Schmid & Henningson (Reference Schmid and Henningson2001)), equation (2.6) can be written entirely in terms of matrix operators as

(2.11a)\begin{equation} (-\mathrm{i} \omega \boldsymbol{\mathsf{M}} + \boldsymbol{\mathsf{L}})\boldsymbol{\mathsf{Q}} = \boldsymbol{\mathsf{B}}_1 \boldsymbol{\mathsf{F}} = \boldsymbol{\mathsf{B}}_1 \boldsymbol{\mathsf{B}}_2 \boldsymbol{\mathsf{R}}, \end{equation}

where the individual matrices are given as

(2.11b)\begin{gather} \renewcommand{\arraystretch}{1.3} \boldsymbol{\mathsf{M}} = \left[\begin{array}{@{}cc@{}} -\Delta & 0\\ 0 & 1 \end{array}\right] ; \renewcommand{\arraystretch}{1.3} \quad \boldsymbol{\mathsf{L}} = \left[\begin{array}{@{}cc@{}} {Re}^{{-}1} \Delta^2 + \mathrm{i} k_x \left(\mathrm{d}^2 \bar{u} - \bar{u} \Delta \right) & 0\\ \mathrm{i} k_z \mathrm{d} \bar{u} & -{Re}^{{-}1} \Delta + \mathrm{i} k_x \bar{u} \end{array}\right]; \end{gather}

(2.11c)\begin{gather} \renewcommand{\arraystretch}{1.3} \boldsymbol{\mathsf{B}}_1 = \mathrm{i} \left[\begin{array}{@{}ccc@{}} k_x \mathrm{d} & -\mathrm{i} k^2 & k_z \mathrm{d} \\ k_z & 0 & - k_x \end{array}\right]; \end{gather}

(2.11d)\begin{gather} \renewcommand{\arraystretch}{1.3} \boldsymbol{\mathsf{B}}_2 = \left[\begin{array}{@{}cccccc@{}} -\mathrm{i} k_x & 0 & 0 & -\mathrm{d} & -\mathrm{i} k_z & 0 \\ 0 & -\mathrm{d} & 0 & -\mathrm{i} k_x & 0 & -\mathrm{i} k_z \\ 0 & 0 & -\mathrm{i} k_z & 0 & -\mathrm{i} k_x & -\mathrm{d} \end{array}\right]; \quad \renewcommand{\arraystretch}{1.3} \boldsymbol{\mathsf{C}} = \dfrac{\mathrm{i}}{k^2}\left[\begin{array}{@{}cc@{}} k_x \mathrm{d} & -k_z \\ - \mathrm{i} k^2 & 0\\ k_z \mathrm{d} & k_x \end{array}\right]; \end{gather}

(2.11e)\begin{gather} \boldsymbol{\mathsf{Q}} = \left[\begin{array}{@{}c@{}} \tilde{V} \\ \tilde{\eta} \end{array}\right];\quad \boldsymbol{\mathsf{R}} =\left[\begin{array}{@{}c@{}} \tilde{R}_{xx} \\ \tilde{R}_{yy} \\ \tilde{R}_{zz} \\ \tilde{R}_{xy} \\ \tilde{R}_{xz} \\ \tilde{R}_{yz} \end{array}\right]; \quad \boldsymbol{\mathsf{U}} = \left[\begin{array}{@{}c@{}} \tilde{U} \\ \tilde{V} \\ \tilde{W} \end{array}\right]; \quad \renewcommand{\arraystretch}{1.3} \boldsymbol{\mathsf{F}} = \left[\begin{array}{@{}c@{}} \tilde{F}_x \\ \tilde{F}_y \\ \tilde{F}_z \end{array}\right]; \end{gather}

where $\boldsymbol{\mathsf{I}}$ is the identity matrix and $\Delta = (\mathrm {d}^2 - k^2)$ is the Laplacian operator, with $\mathrm {d} \equiv \mathrm {d}/\mathrm {d} y$ and $k^2 = k_x^2 + k_z^2$. The velocity–vorticity vector, $\boldsymbol{\mathsf{Q}}$, is related to the Cartesian velocity components, $\boldsymbol{\mathsf{U}}$ as $\boldsymbol{\mathsf{U}}=\boldsymbol{\mathsf{C}}\boldsymbol{\mathsf{Q}}$, where $\tilde {\eta }$ is the wall-normal vorticity component.

Isolating the velocity modes, $\boldsymbol{\mathsf{Q}}$ in (2.11) and rewriting them in Cartesian form, $\boldsymbol{\mathsf{U}}$, yields

(2.12)\begin{equation} \underbrace{\boldsymbol{\mathsf{C}}\boldsymbol{\mathsf{Q}}}_{\boldsymbol{\mathsf{U}}} = \underbrace{\boldsymbol{\mathsf{C}} (-\mathrm{i} \omega + \boldsymbol{\mathsf{M}}^{{-}1} \boldsymbol{\mathsf{L}})^{{-}1} (\boldsymbol{\mathsf{C}}^H \boldsymbol{\mathsf{C}})^{{-}1} \boldsymbol{\mathsf{C}}^H}_{\mathcal{H}(\boldsymbol{k})} \underbrace{\boldsymbol{\mathsf{B}}_2 \boldsymbol{\mathsf{R}}}_{\boldsymbol{\mathsf{F}}}, \end{equation}

where superscript $H$ denotes the Hermitian transpose. Thus we can define a transfer function $\mathcal {H}(\boldsymbol {k})$ to relate the nonlinear Reynolds mode input, $\boldsymbol{\mathsf{F}}$, to the velocity output modes, $\boldsymbol{\mathsf{U}}$. This formulation is identical to the traditional resolvent as discussed in Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013), except for a sign reversal in their equation (2.7b), in element $(2.1)$ of their forcing operator $C^{\dagger}$ (corresponding to our $\boldsymbol{\mathsf{B}}_1$).

The dynamics of the Reynolds stress in (2.7a) can also be written in terms of the substituted normal modes and then rewritten in matrix form for each $(k_x,k_z,\omega )$ as

(2.13a)\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{\mathsf{R}}+\boldsymbol{\mathsf{J}}\boldsymbol{\mathsf{U}} = \boldsymbol{\mathsf{G}}, \end{equation}

where the operator matrices are defined as

(2.13b)\begin{equation} \renewcommand{\arraystretch}{1.3} \boldsymbol{\mathsf{A}} = \mathrm{i}(-\omega + k_x \bar{u}) \left[\begin{array}{@{}cccccc@{}} 1 & 0 & 0 & 2\gamma & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & \gamma & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & \gamma \\ 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right] - {Re}^{{-}1}(\mathrm{d}^2-k^2)\boldsymbol{\mathsf{I}}; \end{equation}

(2.13c)\begin{equation} \renewcommand{\arraystretch}{1.3} \boldsymbol{\mathsf{J}} = \left[\begin{array}{@{}ccc@{}} 2\bar{r}_{xx}\mathrm{i}k_x + 2\bar{r}_{xy}\mathrm{d} & \bar{r}_{xx,y} & 0 \\ 0 & \bar{r}_{yy,y} + 2 \bar{r}_{xy} \mathrm{i} k_x + 2 \bar{r}_{yy} \mathrm{d} & 0 \\ 0 & \bar{r}_{zz,y} & 2\bar{r}_{zz}\mathrm{i}k_z \\ \bar{r}_{xy} \mathrm{i} k_x + \bar{r}_{yy}\mathrm{d} & \bar{r}_{xy,y} + \bar{r}_{xx} \mathrm{i} k_x + \bar{r}_{xy} \mathrm{d} & 0 \\ \bar{r}_{zz} \mathrm{i} k_z & 0 & \bar{r}_{xx} \mathrm{i} k_x + \bar{r}_{xy} \mathrm{d} \\ 0 & \bar{r}_{zz} \mathrm{i} k_z & \bar{r}_{xy} \mathrm{i}k_x + \bar{r}_{yy} \mathrm{d}, \end{array}\right] \end{equation}

with $\gamma =(\mathrm {d}\bar {u}/\mathrm {d} y)/[\mathrm {i}(-\omega + k_x\bar {u})]$; $\overline {(\;)}_{,y} \equiv \mathrm {d}\overline {(\;)}/\mathrm {d} y$, and $\boldsymbol{\mathsf{G}}$ corresponds to the forcing term $\tilde {g}_{ij}$ in (2.7).

The denominator of $\gamma$ expresses the difference between the streamwise phase speed of the normal mode, $c = \omega /k_x$, and the local mean velocity in the channel, $\bar {u}$, and thus $\boldsymbol{\mathsf{A}}$ exhibits a inviscid singularity at the wall-normal location, $y_c$, when the two speeds are matched, such that $\bar {u}(y_c)=c$. This singularity also appears in the large-scale dynamics, within the resolvent operator $\mathcal {H}(\boldsymbol {k})$. In both cases, viscosity (expressed through ${Re}$) resolves this singularity via the generation of a critical layer in the neighbourhood of $y_c$, near where the amplitude of the related ‘critical’ mode reaches its maximum.

2.4. Most amplified velocity modes

In order to obtain the modes shapes for the large- and small-scale structures represented by $\boldsymbol{\mathsf{U}}$ and $\boldsymbol{\mathsf{R}}$, we apply the resolvent framework to the analysis of (2.11) and (2.13). Assuming that the fluctuating Reynolds stresses, $\boldsymbol{\mathsf{R}}$ (and thus the related nonlinear forcing, $\boldsymbol{\mathsf{F}}$) is unknown, the most amplified velocity modes can be modelled by a SVD of the operator, $\mathcal {H}$, linking $\boldsymbol{\mathsf{U}}$ and $\boldsymbol{\mathsf{F}}$, under a kinetic energy norm, following Reddy, Schmid & Henningson (Reference Reddy, Schmid and Henningson1993).

The SVD of the resolvent operator evaluated at a wavenumber triplet $\boldsymbol {k_f}$ can be expressed as the sum of left, $\psi _j$, and right, $\phi _j$, singular vectors, each associated with an ordered sequence of singular values, $\sigma _j$

(2.14)\begin{equation} {\mathcal{H}} = \sum_{j=1}^{\infty} \psi_j \sigma_j \phi_j^H. \end{equation}

The orthonormal singular vectors represent the shape (phase and amplitude) of the velocity (output) and forcing (input) modes that most strongly amplify the kinetic energy associated with the resolvent operator, as given (in vector form) by

(2.15)$$\begin{gather} \tilde{F}_i = \sum_{j=1}^{\infty} \phi_{i, j} \bar{\chi}_{j}, \end{gather}$$

(2.16)$$\begin{gather}\tilde{U}_i = \sum_{j=1}^{\infty} \sigma_{j} \psi_{i, j} \bar{\chi}_{j} = \sum_{j=1}^{\infty} \psi_{i, j} \chi_{j} , \end{gather}$$

where $\sigma _j$ represents the magnitude of the linear amplification of mode $j$, and $\bar {\chi}_{j}$ represents the a priori unknown projection coefficient of the true forcing onto the linearly amplified forcing mode, such that the inputs and outputs of the system are self-consistent. For broadband forcing, $\bar {\chi}_{j}$ is assumed to be unity (i.e. the forcing and principal forcing directions are aligned). Combining the projection coefficient with the linear amplification captured by the singular value, McKeon (Reference McKeon2017) defined an overall ‘velocity weighting factor’, $\chi_{j}= \sigma _{j} \bar {\chi}_{j}$, which conveniently captures the energy of individual disturbance modes, as described below in § 4.

The rank-1 mode of the velocity response modes (i.e. just the first singular mode, $j=1$) has been shown to successfully capture the near-wall cycle and other energetically significant features of wall-bounded flows due to the low rank of the resolvent operator – see McKeon (Reference McKeon2017) for a review. Thus, for the subsequent scale-interaction analysis, the mode shapes of $\tilde {U}_i$ will be modelled by only $\psi _{i,1}$.

In order to obtain the Reynolds modes, $\tilde {R}_{ij}$ corresponding to the velocity modes $\tilde {U}_i$, there are two possible approaches. In Jacobi et al. (Reference Jacobi, Chung, Duvvuri and McKeon2021), they assumed that the nonlinear forcing to the Reynolds stress, $\boldsymbol{\mathsf{G}}$, appearing in (2.13) was uncorrelated with the scale interactions (a variation on the RNL model)

(2.17)\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{\mathsf{R}}+\boldsymbol{\mathsf{J}}\boldsymbol{\mathsf{U}} \approx 0 \end{equation}

and thus expressed $\boldsymbol{\mathsf{R}}$ in terms of the $\boldsymbol{\mathsf{U}}$ modes obtained above via the resolvent analysis, according to

(2.18)\begin{equation} \boldsymbol{\mathsf{R}} \approx{-}(\boldsymbol{\mathsf{A}}^{{-}1} \boldsymbol{\mathsf{J}})\boldsymbol{\mathsf{U}}. \end{equation}

However, in the previous analysis, they did not actually solve for the resolvent modes via the SVD. Rather, they approximated the $\tilde {U}_i$ and $\tilde {R}_{ij}$ mode shapes analytically, assuming they exhibited the self-similar resolvent form identified in Dawson & McKeon (Reference Dawson and McKeon2019). In the present analysis, we adopted this same general procedure, but we represented the $\boldsymbol{\mathsf{U}}$ mode shape with the first singular mode, $\psi _{i,1}$, obtained from the numerical solution of (2.11), without any additional analytical modelling.

Upon substitution of the rank-1, large-scale mode shapes, $\psi _{i,1}$, into (2.18), we denote the resulting rank-1 Reynolds stress mode shapes as $\xi _{i,1}$. Note that these Reynolds stress modes are not orthonormal.

An alternative approach to obtaining the $\boldsymbol{\mathsf{U}}$ and $\boldsymbol{\mathsf{R}}$ mode shapes, in which $\boldsymbol{\mathsf{G}}$ is not neglected and instead the resolvent framework is applied directly to the $\boldsymbol{\mathsf{R}}$ modes, is discussed in Appendix B.

3.1. Phase-difference calculations

The average phase differences between the Fourier modes, shown in figure 2(c), were obtained directly from the ensemble-averaged cross-spectrum of the DNS, $\hat {u}^* \widehat {u^2}$. The phase differences were then unwrapped to obtain a continuous phase-difference profile. However, it is important to note that the phase difference between the two Fourier modes cannot be obtained by subtracting the phases of the separate, ensemble-averaged mode shapes in figure 2(d,e) due to the non-zero covariance between the velocity and stress signals, i.e. $\langle \hat {u}^* \widehat {u^2} \rangle = \langle \hat {u}^* \rangle \langle \widehat {u^2} \rangle + \text {cov}(\hat {u}^*,\widehat {u^2})$, where $\langle {\cdot } \rangle$ represents ensemble averaging. This covariance problem applies more generally to any mode eduction by ensemble averaging, including SPOD. The phase difference observed between individual, ensemble-averaged modes will not correspond to the phase difference implicit in the amplitude modulation coefficient (2.1) that we are seeking to predict. For this same reason, the shape of the individual modes predicted by the resolvent cannot be compared with ensemble-averaged modes educed from measurements or computations, since those educed modes do not represent the ensemble phase difference, whereas the transfer function between the resolvent modes does.

As expected from the phase interpretation of the amplitude modulation coefficient, the phase profile extracted from the cross-spectrum starts at $\Delta \phi = 0$ at the wall, where the velocity and stress are in phase, and decreases monotonically to $-{\rm \pi}$ far from the wall, where the velocity and stress are out of phase. We then compared this phase-difference profile with the difference predicted by the resolvent mode.

The phase difference, $\Delta \phi = \phi _{\xi _{1,1}}- \phi _{\psi _{1,1}}$ between the large-scale velocity and stress modes, calculated via the resolvent framework, was obtained by unwrapping the phases of the $\psi _{1,1}$ and $\xi _{1,1}$ modes from the wall outwards before subtraction and is illustrated in figure 2(c) in the solid line.

The phase calculation from the DNS cross-spectrum and the analysis of Jacobi et al. (Reference Jacobi, Chung, Duvvuri and McKeon2021) both indicated that the phase difference, $\Delta \phi$, is negative across the wall layer, consistent with many other empirical studies. By contrast, the resolvent phase difference calculated here is positive very near the wall and then changes sign farther out in the wall layer. At the location of the outer spectral peak, $y_{op}$, previous studies and the DNS results indicated a phase difference of around $-{\rm \pi} /2$, consistent with the zero crossing of the amplitude modulation coefficient, $R$ (since the amplitude modulation coefficient is related to the phase difference, heuristically speaking, via $R \lesssim \cos {(\Delta \phi )}$ as noted above; see in § 2.1 for details), but the resolvent modes predicted a positive phase of approximately $+{\rm \pi} /4$.

The discrepancy between the resolvent and DNS phase calculations originates immediately at the wall, where previous studies reported a negligible phase difference and the single phase-speed resolvent predicted a positive difference. This phase disagreement at the wall is likely a result of the differences between the Fourier and resolvent mode shapes in this region, the two most prominent of which are the downstream inclination angles and near-wall modal amplitudes, as visible in figure 2. However, the discrepancy in inclination angle itself is not necessarily the primary cause for concern, since both the velocity and stress modes are similarly inclined and we are interested only in the relative phase difference between the two. The key concern is actually the discrepancy in amplitude. The resolvent mode exhibits a maximum amplitude at the location of the VLSM peak, where its phase speed matches the local mean velocity (the critical layer), and thus its amplitude closer to the wall is much lower than that of the corresponding spatial Fourier mode at the same wavenumber. This means that, from the perspective of a modal decomposition of the flow, the particular resolvent mode centred at the VLSM peak, $\tilde {U}(y,k_x,k_z,c_{op})$, is a poor analogue for the spatial Fourier mode educed from the DNS data, $\hat {u}(y,k_x)$.

3.2. The problem with a single resolvent mode model

The fundamental problem with comparing a phase model based on a single resolvent mode with the Fourier DNS modes is that the spatial Fourier mode is really an ensemble average of the full, spatio-temporal Fourier mode, $\hat {u}(y,k_x,k_z,\omega )$, over frequency (or, equivalently, phase speed) and spanwise wavenumber. Therefore, the spatial Fourier mode incorporates energy from a variety of phase speeds at a given wall-normal location. The range of different convection velocities, $c$, associated with a single wavenumber $k_x$ of spatial Fourier modes can be observed from the non-negligible width of two-dimensional, temporal/spatial energy spectra (Lehew, Guala & McKeon Reference Lehew, Guala and McKeon2011). Thus, we can think of the relationship between the spatial Fourier mode, $\hat {u}(y,k_x)$, and the resolvent modes, $\tilde {U}_1(y,k_x,k_z,c)$, in terms of an ensemble average

(3.1)\begin{equation} \hat{u}(y,k_x) \approx k_x \iint _{c k_z} \varGamma(k_z,c) \tilde{U}_1(y,k_x,k_z,c) \, {\rm d} k_z \, {\rm d}c, \end{equation}

where $\varGamma (k_z,c)$ is a weighting function, representing the (unknown) joint probability density function (p.d.f.) of phase speeds and spanwise wavenumbers of the resolvent modes at a fixed wall-normal location. By utilizing a single phase speed for our resolvent mode, above, we assumed that (i) a single resolvent mode with phase speed, $c_{op}$, dominated the modes of all other phase speeds; and (ii) that this dominance applied uniformly across all $y$-locations. However, these assumptions depend on the shape and amplitude of the resolvent modes and are not satisfied.

McKeon, Sharma & Jacobi (Reference McKeon, Sharma and Jacobi2013) identified two distinct regimes of resolvent-generated modes: attached modes, which can appear as either wall modes or critical modes, and detached critical modes. Critical modes are distinguished from wall modes by expressing a localized amplitude maximum at the wall location, $y_c$, close to where the local mean velocity matches the phase speed, $\bar {u}(y_c) =c$, known as the critical point. Attached modes are distinguished from detached by expressing a high amplitude even very near the wall, thus exemplifying Townsend's concept of exerting a ‘footprint’ on the wall. Because the detached modes are localized away from the wall, they obtain a self-similar universal form for a range of phase speeds $16 < c < U_0 - 6.15$ (Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013) that allows for the asymptotic modelling of the mode shape noted above. But this also means that detached, self-similar modes have very low amplitude near the wall, as in the case of our single critical mode considered above with a critical layer around $y_{op}$.

Figure 3(a) shows the amplitude map of the principal streamwise velocity mode, $\psi _{1,1}(y)$, with varying convective velocity, $c$, for the streamwise VLSM structures calculated by the resolvent approach. For high convective velocities, we observe the detached, self-similar critical modes that were utilized in the analytical modelling, each centred on the wall-normal location matching its local mean velocity (marked by the black solid line). At lower convective velocities, the critical modes begin to attach to the wall, and lower still eventually lose their ‘critical’ characteristic shape and appear as attached wall modes.

Figure 3. (a) Amplitude map of the VLSM mode, $\psi _{1,1}$, across the half-channel as a function of phase speed, $c$. The blue dashed line marks the mean velocity (and thus critical phase speed) at the outer spectral peak, $y_{op}$, and thus traces along the amplitude profile of the mode with $c_{op} = 0.71 U_0$. Due to the very low amplitude near the wall, it is clear that this mode is detached. The other coloured lines indicate large-scale modes with different phase speeds, ranging from detached to attached. These 4 modes are illustrated in the cartoon in (b). The lowest phase speed (grey), representing the mean velocity at around $y^+ \approx 10$, generates an attached wall mode. As the phase speed increases, the resolvent modes become critical and eventually detach.

Four sample modes taken from a near-wall location in figure 3(a) are sketched in figure 3(b) to illustrate how they disconnect from the wall region with increasing convective velocity. The detached (blue) mode is the mode used in the single velocity analysis above; it convects at the local mean velocity at the location of $y_{op}$. Because it is a critical mode, its amplitude will dominate the amplitude of any other resolvent mode at the same wall-normal location. But it will not dominate over other resolvent modes closer to the wall. In fact, nearer to the wall, the amplitude of the resolvent mode defined by phase speed $c_{op}$ is negligible compared with other phase speeds. Thus the single mode centred on $y_{op}$ will not describe the statistically representative (or ‘average’) mode shape in the region close to the wall.

The fact that the critical resolvent mode centred at the outer peak is not, in fact, the dominant mode nearer to the wall creates a significant problem for computing a meaningful phase difference profile. When the phase difference from this outer-peak-centred resolvent mode was unwrapped, it incorporated phase information from locations where it was not dominant, and thus where its associated phase was also not representative. But phase unwrapping depends on the phase values everywhere along the unwrapping path in order to be reliable. If the phase near the wall is based on a mode that does not actually represent the flow, then the unwrapped phase farther from the wall where that mode is indeed representative will be contaminated. For the unwrapped phase to be physically meaningful, it must reflect the contributions of only significant, representative modes along the entire unwrapping path, from the wall to any point of interest. (This is a problem specific to the need to unwrap phase profiles and is not relevant for standard superposition analysis. The analytical study of Jacobi et al. (Reference Jacobi, Chung, Duvvuri and McKeon2021) avoided this problem by assuming a zero phase difference at the wall and a monotonic unwrapping out to the neighbourhood of the critical layer at $y_{op}$, without verification.)

In order to predict a meaningful profile of the phase difference, we must incorporate the phase information from different resolvent modes, each one dominant and representative of the flow at a different wall-normal height, instead of just a single mode, in order to properly unwrap the phase. The challenge is therefore to determine how to construct a statistically representative composite mode.

4.1. Cumulative phase composition

In order to calculate the phase difference between different modes at a particular wall-normal location, $y$, the phase of each mode must first be unwrapped at all $y$-positions, from the wall out to that location. Without unwrapping, the phase difference is meaningless, e.g. a wrapped phase difference of $+{\rm \pi} /2$ might actually be associated with a true phase difference of $-3 {\rm \pi}/2$ that wrapped around to the positive branch of the argument function. In the wrapped state, it is impossible to establish the lead or lag of different modes. Therefore, the correct value of the phase difference is entirely dependent on the validity of its unwrapping from the wall.

This presents a problem if we consider a piecewise composition of spatially localized mode shapes (like resolvent modes or SPOD modes), because the individual (constituent) modes used for the pieces of the composition cannot be unwrapped separately and then combined. Rather, the phase of the overall composite mode must be unwrapped as part of the composition process, itself.

In order to accomplish this simultaneous mode composition and unwrapping, we consider piecing together a set of dominant resolvent modes, each with phase speed $c^*$ resulting in a maximum (weighted) amplitude at location $y^*$, such that the overall composite mode includes phase and amplitude information only from the neighbourhoods of these amplitude maxima, as sketched in simplified form in figure 4(a–c). In other words, one resolvent mode at each $y^*$ location (defined by its phase speed $c^*$) is assumed to contribute predominantly to the composite mode at that same location. Therefore, holding wavenumbers constant, we assume that $(y^*,c^*(y^*))$ identifies the phase speed that maximizes the weighted mode amplitude at any given $y^*$ location, shown in figure 4(a). And, for each constituent velocity or stress mode with phase speed $c^*$, the phase profile is given by $\phi (y,k_x,k_z,c^*)$, shown in figure 4(b).

Figure 4. A heuristic illustration of the phase composition method described in (4.1) for the $\psi _{1,1}$ mode shape. At each height $y_i^*$, we identify the energetically dominant resolvent mode with phase speed $c_i^*$ from among all the different phase speeds. The amplitude profiles for three such dominant modes are shown in (a). The phase profiles for each of these modes are unwrapped from the wall, individually, in (b). The derivative of each phase profile, in the neighbourhood of its maximum amplitude, is then composited piecewise and integrated from the wall, in order to obtain an unwrapped phase profile that preserves the phase information from the dominant modal contributions at each height, as shown in (c).

Now, we can define the composite phase profile, $\langle \phi (y,k_x,k_z) \rangle$, for a piecewise composition of modes with different phase speeds, $c^*$ at each wall-normal location, $y^*$, according to

(4.1)\begin{equation} \langle \phi(y,k_x,k_z) \rangle = \phi(0,k_x,k_z,c^*_0)+\int_{0}^{y} \frac{\partial}{\partial y} \phi(y,k_x,k_z,c^*) \big|_{y=y^*} \, {{\rm d} y}^*, \end{equation}

where $c^*_0 \equiv c^*(0)$ is the phase speed corresponding to $y^*=0$, at the wall. This anti-derivative accumulates the phase of the different constituent modes as they are integrated together to form the composite mode, such that the phase of the resulting composite mode varies continuously and captures the slope of each constituent piece, illustrated in figure 4(c). The derivative in the integrand can be approximated numerically in the neighbourhood of the peak amplitude location, $y^*$, to any order of accuracy by writing a series expansion of $\phi (y,k_x,k_z,c^*)$ near $y^*$. Thus, we can use (4.1) to calculate a meaningful phase profile of a composite mode that is constructed from the dominant resolvent modes at all (or a subset of all) wall-normal positions.

Now, the problem is how to determine which are the dominant resolvent modes, or, in other words, what is the phase speed, $c^*$, of the dominant mode at each wall-normal location.

4.2. Determining the dominant resolvent modes

To compare the relative importance of modes of different convective velocities, $c$, at a given $y$ location, we consider two weighting approaches based on kinetic energy. Following the definitions of the velocity modes in (2.16), the streamwise spectral energy density of the rank-1 approximated VLSM modes, $\varPhi _1$, can be written in two different forms, depending on assumptions about the forcing

(4.2)\begin{equation} \varPhi_1 = \tilde{U}_1^* \tilde{U}_1 = \begin{cases} \sigma_{1}^2 \, |\psi_{1, 1}|^2, & \text{if}\ |\bar{\chi}_{1}|^2=1, \quad \text{broadband forcing}, \\ |\chi_{1}|^2 \, |\psi_{1, 1}|^2, & \text{unknown forcing}. \end{cases} \end{equation}

Assuming broadband forcing, the forcing projection term, $\bar {\chi}_1$ can be set to unity, and the upper expression depends on the singular value, $\sigma _1^2$, which expresses the linear amplification of the mode, and on $\psi _{1,1}$ which expresses how that amplification is distributed over the mode's amplitude profile. This version of the kinetic energy distribution over the mode shape has the advantage of being entirely predictable from within the resolvent itself, via the singular value.

Relaxing the broadband forcing assumption, the lower expression utilizes a general velocity weighting, $\chi_1^2$, to express the total (linear and nonlinear) weighting necessary to render the output mode shapes consistent with the unknown, nonlinear forcing. The nonlinear forcing effect was also pointed out in McKeon (Reference McKeon2017) and the idea of the colour of forcing is explored in detail in Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021). This weighting is not, a priori, known and must be modelled, but would, in principle, have the advantage of greater fidelity to true turbulent flows. We consider both approaches, in order to identify the phase speed, $c^*$, of the dominant resolvent mode.

4.3. Resolvent weighting: linear amplification

We adopt the broadband forcing assumption first, for simplicity. Figure 5(a) shows the map of the spectral energy density, $\varPhi _1$, of the VLSM modes under this assumption, as a function of $y$ and $c$. The particular convection velocity, $c^*$, of the VLSM that contributes most strongly to the energy of the system at height $y$, is then given by tracing along the peak of this map, and can be written formally as

(4.3)\begin{equation} c^*(y^*,k_x,k_z) = \underset{c}{\arg\max} [\sigma_1^2(y^*,k_x,k_z,c) |\psi_{1,1}(y^*,k_x,k_z,c)|^2], \end{equation}

which is illustrated by the solid red line over the map. This function $c^*(y)$ indicates the equivalent convective velocity profile for a representative VLSM mode that can be constructed piecewise from the dominant mode contributions at each height, $y$, with each piece convecting at a different velocity. The resulting composite modes for the VLSM are shown in figure 5(b); the corresponding small-scale (Reynolds stress) mode shapes, $\xi _{1,1}$, calculated from the same piecewise procedure as the large scales, are shown in figure 5(c).

Figure 5. (a) Spectral energy density profiles, $\sigma _1^2 |\psi _{1, 1}|^2$, under the assumption of broadband forcing given in (4.2). The red line represents the convection velocity, $c^*$, defined in (4.3) used to construct the composite modes; (b) the composite velocity mode, $\psi _{1, 1}$; and (c) the composite Reynolds stress mode $\xi _{1, 1}$ defined by the convective velocity profile in (a); the normalization and phase line markings are the same as figure 2. (d) The phase difference, $\Delta \phi$, calculated from the composite modes is shown by the black line, contrasted with the phase from the DNS.

The piecewise modes can be compared with the single speed mode shown above in figure 2(a,b). The Reynolds modes in both cases show an initial upstream inclination and then, farther from the wall, a downstream inclination, leading to the same, problematic, positive phase difference between the large and small scales at the wall.

The profile of the phase difference between the large and small scales from the piecewise modes is shown in figure 5(d). Although the sense of the phase difference is still incorrect very near the wall, the trend appears slightly more realistic than the single mode calculation, at least farther from the wall. Nevertheless, weighting based on the linear amplification predicted by $\sigma _1^2$ still does not seem sufficient to capture a representative signature of the VLSMs, and thus we turn to the second weighting option, involving the general velocity weighting factor, $\chi_1$.

4.4. Resolvent weighting: modelling nonlinear weights

Because the general velocity weighting cannot be determined endogenously from the resolvent calculation itself, it requires a model. A number of previous studies have considered how to model this weighting in order to align the spectral energy density predicted by the rank-1, critical mode resolvent approximation, $\varPhi _1(y,k_x,k_z,c)$, with the true spectral energy density, $\varPhi (y,k_x,k_z,\omega )$. Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013) introduced a phase-speed-dependent weighting function, $W({c})$, and constructed an optimization problem to correct the integrated form of $\varPhi _1$, whereas Moarref et al. (Reference Moarref, Jovanović, Tropp, Sharma and McKeon2014) applied a similar optimization approach to a matrix constructed from the velocity weighting factors, $\chi_j$, directly.

Instead of utilizing an optimization approach, in this analysis, we will adopt an empirical/heuristic argument to model the velocity weighting factors, $\chi_j$, in order to properly weight the contributions of VLSM resolvent modes at different phase speeds for a given wall-normal location.

We begin with the empirical observation (Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013) that the pre-multiplied, streamwise spectral energy density of wall-bounded turbulence, integrated over $\omega$, and evaluated at the $y$ location of the critical layer, $y_{c'}$, is similar in form to the low rankness measure of the rank-1 singular modes, $S(k_x,k_z,c') \equiv \sigma _1^2/\sum _{j=1}^{\infty } \sigma _j^2$, evaluated at the corresponding convection velocity, $c'$ for that same critical layer

(4.4)\begin{equation} k_x k_z \int_\omega \varPhi(y_{c'},k_x,k_z,\omega) \, {\rm d}\omega \sim S(k_x,k_z,c'), \end{equation}

where the left-hand side is a function of $y_{c'}$ which is the specific location compatible with the specific phase speed $c'$ on the right-hand side, $\bar {u}(y_{c'}) = c'$. We assume that $(k_x,k_z,c')$ are fixed on both sides of this similarity expression.

The spectral energy density is then expressed by its rank-1 approximation including the velocity weighting

(4.5)\begin{equation} \varPhi_1(y,k_x,k_z,c) = \tilde{U}_1 \tilde{U}_1^* = |\chi_{1}(k_x,k_z,c)|^2 |\psi_{1,1}(y,k_x,k_z,c)|^2, \end{equation}

and this is substituted into the similarity statement, changing the integration from frequency to phase speed, $d\omega = k_x \, {\rm d}c$, to obtain

(4.6)\begin{equation} k_x^2 k_z \int_c |\chi_{1}(k_x,k_z,c)|^2 |\psi_{1,1}(y_{c'},k_x,k_z,c)|^2 \, {\rm d}c \sim S(k_x,k_z,c'), \end{equation}

where we note that the velocity weighting factor $|\chi_{1}(k_x,k_z,c)|^2$ appears similar to the $W({c})$ weighting function constructed in Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013), with the important difference that the current factor is wavenumber dependent.

Because the mode shape dependence of the kinetic energy is already accounted for via $\psi _{1,1}$, we focus on just the convection velocity effect in the weighting term. Therefore, we integrate the similarity statement over wall-normal position, $y_{c'}$

(4.7)\begin{equation} k_x^2 k_z \int_c |\chi_{1, 1}(k_x,k_z,c)|^2 \left[ \int_{y_{c'}}|\psi_{1}(y_{c'},k_x,k_z,c)|^2 \, {\rm d}y_{c'} \right] \, {\rm d}c \sim \int_{y_{c'}} S(k_x,k_z,c') \, {\rm d}y_{c'}, \end{equation}

and exploit the unit normalization of the integrated singular modes to eliminate the bracketed term on the left-hand side

(4.8)\begin{equation} k_x^2 k_z \int_c |\chi_{1}(k_x,k_z,c)|^2 \, {\rm d}c \sim \int_{y_{c'}} S(k_x,k_z,c') \, {\rm d}y_{c'}. \end{equation}

To simplify the right-hand side, we change the wall-normal integral into an integral over convection velocity, according to ${\rm d}y_{c'} = ({{\rm d}\bar {u}}/{{\rm d} y})^{-1}|_{y_{c'}} \,{\rm d}c'$ and obtain

(4.9)\begin{equation} k_x^2 k_z \int_c |\chi_{1}(k_x,k_z,c)|^2 \, {\rm d}c \sim \int_{c'} S(k_x,k_z,c') \, \left(\frac{{\rm d}\bar{u}}{{\rm d} y} \right)^{{-}1}_{y_{c'}} \,{\rm d}c'. \end{equation}

Note here that $({{\rm d}\bar {u}}/{{\rm d} y} )^{-1}_{y_{c'}} = \left ({{\rm d} y}/{{\rm d}\bar {u}} \right )_{c'}$ and thus is a function of $c'$. The integrals on both sides of the similarity can be written over $c$ and thus the integrands themselves must be similar, such that

(4.10)\begin{equation} |\chi_{1}(k_x,k_z,c)|^2 \sim \frac{1}{k_x^2 k_z} \left(\frac{{\rm d} y}{{\rm d}\bar{u}} \right)_{c} S(k_x,k_z,c), \end{equation}

and we obtain an expression for the otherwise unknown velocity weighting. (The proportionality factor is unimportant due to the ${\arg \max }$ operation.)

Based on this empirical/heuristic argument, when seeking the particular convection velocity $c^*$ that corresponds to the VLSM with the largest energetic contribution at a particular $y$ location, we solve

(4.11)\begin{align} c^*(y,k_x,k_z) &= \underset{c}{\arg\max} [ |\chi_1(k_x,k_z,c)|^2 |\psi_{1,1}(y^*,k_x,k_z,c)|^2 ] \end{align}

(4.12)\begin{align} & \approx \underset{c}{\arg\max} \left[ \frac{1}{k_x^2 k_z} \left(\frac{{\rm d} y}{{\rm d}\bar{u}} \right)_{\bar{u}=c} \frac{\sigma_1^2}{\sum_{j=1}^{\infty} \sigma_j^2} |\psi_{1,1}(y,k_x,k_z,c)|^2 \right]. \end{align}

Figure 6(a) shows the map of the energy distribution of VLSM modes with this weighting model, as a function of $y$ and $c$, where the magenta solid line indicates the dominant convective velocities, $c^*(y)$. Note that the convective velocity of the VLSM approaching the wall is predicted to be $0.59 U_0$, precisely the value that would eliminate the problematic phase difference at the wall observed in § 3.

Figure 6. Caption entries follow from figure 5 with (a) the spectral energy density profile map and optimal convection velocity defined by (4.11); (b,c) the constructed piecewise mode shapes for $\psi _{1, 1}$ and $\xi _{1,1}$, respectively; and (d) the resulting phase difference profile between scales, compared with DNS results.

The large- and small-scale mode shapes generated by the piecewise profile of convective velocities are shown in figure 6(b,c), respectively. The differences between these modes and previous mode shapes are quite significant. In particular, the Reynolds mode no longer exhibits a change in inclination angle near the wall, and as a consequence, the weighting model results in a negligible phase difference between large and small scales at the wall itself, consistent with experimental observations and with the assumptions of the previous analytical approach.

The profile of the phase difference, $\Delta \phi$, is shown in figure 6(d) and compares reasonably well with that calculated from the DNS channel data at matched Reynolds number. Unlike the single mode or linear-weighted, piecewise mode approaches, the sign of the phase difference appears accurate across the entire near-wall region. Using the new $\chi_j$ weighting therefore results in a reasonably close estimate of the inter-scale phase profile.

However, the agreement between the Fourier and composite modes is still not perfect. The inclination of the composite velocity mode is still much steeper than the corresponding Fourier spatial mode. And near the wall, the composite mode seems to under-predict the magnitude of the phase lag, while over-predicting it farther from the wall. This discrepancy is likely a result of the single-wavenumber simplification at the foundation of the modelling approach, in which we inferred from the cross-spectral power in figure 1(c) that a narrow range of wavenumbers dominates the inter-scale coupling, and thus we excluded terms with wavenumbers beyond the specific scale of interest from (2.6a). But this simplification does not apply near the wall, where the range of scales participating in scale interactions broadens significantly. Moreover, we assumed that only a single spanwise wavenumber, $k_z$, dominated the flow. It is also important to note that, because the weighting formula was derived from the streamwise energy spectrum, its application to other velocity components is less certain, and other weighting models may be required, as discussed in Appendix D.

Nevertheless, the weighted profile reproduces the wall-normal location where the phase difference is $-{\rm \pi} /2$ to within 20 %, which is less than the $40\,\%$ relative variation in zero-crossing location due to varying the filter cutoff reported by Mathis et al. (Reference Mathis, Hutchins and Marusic2009). The $-{\rm \pi} /2$ location of the predicted phase profile is important because the resolvent model prediction of the phase difference profile, $\Delta \phi (k_x)$ represents a wavenumber-specific contribution to the amplitude modulation coefficient, $R$, which averages phase-difference contributions over all wavenumbers, $R = \langle \cos {(\Delta \phi (k_x))} \rangle$. Thus, the $-{\rm \pi} /2$ phase location corresponds to the zero-crossing point of the amplitude modulation coefficient and the VLSM energetic peak at around $y^+= 3.9Re^{1/2}$. And the profile of the amplitude modulation coefficient, in turn, plays an important role in efforts at inner-outer wall modelling. In particular, Mathis, Hutchins & Marusic (Reference Mathis, Hutchins and Marusic2011a) showed that their modelling coefficient, $\beta$, which regulates the magnitude of influence of the outer large-scale fluctuations, is basically identical to $R$. Therefore, the resolvent model could be used to generate the $\beta$ coefficient for inner–outer models, independent of experimental data, after which it can be incorporated into large-eddy simulation (LES) wall models that use the inner–outer modelling formulation, like in the work of Inoue et al. (Reference Inoue, Mathis, Marusic and Pullin2012). Or, the inner–outer model could be rewritten in the spectral domain as presented in Baars, Hutchins & Marusic (Reference Baars, Hutchins and Marusic2016b), on a wavenumber-by-wavenumber basis, in which case the wavenumber-specific phase difference could be employed directly.

Because of the importance of identifying the $-{\rm \pi} /2$ phase crossing location, we explored the effect of different choices for the dominant wavenumbers, $k_x$ and $k_z$, on this point of the phase profile.

4.5. Wavenumber and Reynolds number dependence

The above analysis showed that a composite VLSM mode based on appropriately weighted convective velocity contributions can reasonably predict the phase-difference profile between large- and small-scale motions. But that analysis studied only a single VLSM scale, with fixed spatial wavenumbers $(k_x, k_z) = (0.75, 6)$ corresponding to the energetic outer peak in the premultiplied energy spectral density map at a single Reynolds number, ${Re} = 5200$. Other large-scale motions also contribute to the modulation between large and small scales, and thus we examine how the phase difference varies with wavenumber, for a fixed Reynolds number, and also how changing the Reynolds number affects the phase profile.

Figure 7(a) shows the phase difference evaluated at the outer spectral peak location, $\Delta \phi (y_{op})$, over a wide range of streamwise and spanwise wavenumbers corresponding to different isolated scales, overlaid with the contour lines of the premultiplied streamwise spectral energy density. The phase difference is negative for all scales, indicating that the small-scale envelope of fluctuations always (spatially) leads the large scale, irrespective of the wavenumbers. However, the magnitude of the phase lead varies with wavenumber; for a band of large scales with $k_x \sim {O}(1)$, the phase lag at the outer spectral peak is approximately $-{\rm \pi} /2$, consistent with experimental findings. But the magnitude of the lead tends to increase with increasing $k_x$, which means that the $-{\rm \pi} /2$ point between the large and small scales moves closer to the wall, in outer units. We visualize this shift in the $-{\rm \pi} /2$ point directly in figure 7(b), which shows the relative error in the predicted location of the $-{\rm \pi} /2$ phase shift vs the outer spectral peak location, $y_{op}$. As noted above, this location corresponds to the $R=0$ point in the amplitude modulation analysis and can vary by as much as $40\,\%$ due to changing filter cutoffs (Mathis et al. Reference Mathis, Hutchins and Marusic2009). Nevertheless, there is a significant band of large-scale motions across all spanwise wavenumbers which predict the crossing point to within $5\,\%$.

Figure 7. (a) The phase difference $\Delta \phi$ of the composite mode at the location of the outer energy peak $y_{op}$ for varying streamwise and spanwise wavenumbers. The black lines indicate the contour levels [0.1, 0.3, 0.6] of the two-dimensional pre-multiplied energy spectrum $k_xk_z\phi _{uu}(k_x, k_z)$ of ${Re}=5200$ DNS channel flow. The red circle marks the wavenumber, $(k_x, k_z) = (0.75, 6)$, used for the composite modes above in § 3 and § 4. (b) The relative error in the predicted location of the $-{\rm \pi} /2$ phase shift compared with the outer energy peak location, $y_{op}$, across all the wavenumbers. The relative error for the composite mode case (red circle) is $19\,\%$, while a slightly smaller LSM, $(k_x, k_z) = (1, 3)$, (blue square) exhibits an error of only $4\,\%$. The predicted mode shapes between these two cases were not significantly different.

Mathis et al. (Reference Mathis, Hutchins and Marusic2009) showed that as the filter-cutoff wavenumber increased, the zero crossing in $R$ moved away from the wall (and the amplitude of the correlation increased generally). Physically, by increasing the filter-cutoff wavenumber, fewer large scales are included in the small-scale signal. This means that the large- and small-scale signals diverge in phase more slowly (moving away from the wall), because the large scales can exert more influence over a signal containing few energetic large-scale signals, resulting in the outward shift of the zero-crossing location. Decreasing the filter-cutoff wavenumber puts more large scales in the small-scale signal, resulting in less influence of the large scales over other large scales, a faster divergence of phase (away from the wall), and an inward shift of the zero crossing. In other words, the trend in the zero crossing can be explained in terms of the greater ease with which large scales can modulate smaller scales compared with modulating other large scales.

In the current model, the predicted trend in zero crossing is different, because the current model only examines a single scale, without any cutoff filter. Thus the current model is more like the filtered case with relatively large scales appearing in both signals, in which case the divergence happens more quickly due to the weaker ability of one signal to modulate the other. Therefore, for a single-wavenumber resolvent mode, increasing the wavenumber means that equivalent-sized smaller scales are interacting with each other, and are less able to influence one another than larger, energetic scales would be, thus resulting in a faster phase divergence (away from the wall) and a $-{\rm \pi} /2$ shift closer to the wall. The rate of phase divergence of the resolvent model therefore offers a measure of the degree of influence between velocity and stress modes that is scale specific, and supports the idea that targeting lower wavenumber pairings of velocity and stress can result in faster and more effective modulation between the two.

The streamwise wavenumber of the large scales, $k_x$, exerts a much stronger influence on the phase profile behaviour than the spanwise wavenumber, $k_z$, although the best predictions for the $-{\rm \pi} /2$ phase crossing occur for the region aligned with $k_x \approx 0.65 k_z^{1/2}$ or $\lambda _x \approx 3.8 \lambda _z^{1/2}$. Jimenez & Hoyas (Reference Jimenez and Hoyas2008) reported that the two-dimensional spectral contours of the streamwise velocity are aligned along $\lambda _x \approx y^{-1} \lambda _z^2$ curves, representing streamwise-oriented, conical velocity structures. Therefore, the vortex structures associated with the best prediction of the phase profile are also conical velocity structures, but wider and squatter than those associated with the energetic peak. The additional width may represent the statistical effect of streamwise meandering of the large scales, which is not captured by the spanwise homogeneous construction of the resolvent operator. The fact that the spanwise and streamwise dimensions scale disproportionately also indicates that the predominant scales involved in the phase profile construction are not geometrically self-similar, which is possibly a low Reynolds number effect as noted by Deshpande et al. (Reference Deshpande, Chandran, Monty and Marusic2020), although this may also indicate limitations to the chosen weighting, which is based on streamwise energy alone.

The effect of the Reynolds number on the phase-difference prediction was also examined, specifically in the context of the $y$-location where the phase difference between scales equal $-{\rm \pi} /2$ which roughly corresponds to the zero-crossing location of the amplitude modulation coefficient, as noted above, and to the location of the outer spectral energy peak, given empirically as $y_{op}^+ \approx 3.9{Re}^{1/2}$.

To see if this empirical trend is obeyed by the predicted phase differences, we calculated the $y$-location of the $-{\rm \pi} /2$ phase difference against Reynolds number, shown by the red circles in figure 8 for the wavenumber pair used in the earlier sections. Although the predicted trend is consistent with an increasing zero-crossing location with Reynolds number, the prediction appears to saturate much more quickly than experiments. This is likely due to the fact that the resolvent includes only a single wavenumber, $k_x$, for the large scales, whereas the true $R(y)$ signal (on which the zero-crossing empirical observation is based) depends on an integral across all such wavenumbers, as shown in § 2.1. Including different (or multiple) streamwise and spanwise wavenumber contributions, perhaps in self-similar triadic hierarchies as demonstrated in Sharma, Moarref & McKeon (Reference Sharma, Moarref and McKeon2017), would introduce additional phase differences into that integral. In particular, higher values of $k_x$ and lower values of $k_z$ tend to be associated with larger phase differences at a fixed $y$-location, as shown in figure 7(a), which means that the $-{\rm \pi} /2$ location would tend to shift closer to the wall, as noted above. Therefore, the exclusion of these smaller scales would tend to result in exaggerating the wall-normal location of the $-{\rm \pi} /2$ phase, thus possibly explaining why the single-wavenumber model trend differs from empirical trends that are based on a correlation coefficient.

Figure 8. Wall-normal height where the phase difference, $\Delta \phi = -{\rm \pi} /2$ vs Reynolds number. The discrete markers denote resolvent predictions from the current model using mean flow quantities from DNS channel flows at: ${Re} =180, 550, 1000, 2000$ from Hoyas & Jiménez (Reference Hoyas and Jiménez2008) and $5200$ from Lee & Moser (Reference Lee and Moser2015); the red circles are calculated for wavenumbers $(k_x, k_z) = (0.75, 6)$ (corresponding to the red circle in figure 7); the blue squares are calculated for $(k_x, k_z) = (1, 3)$ (the blue square in figure 7). The black dashed line represents the relation of $3.9{Re}^{1/2}$ from Mathis et al. (Reference Mathis, Hutchins and Marusic2009), and red dashed line is $0.42{Re}^{3/4}$ inferred from Klewicki et al. (Reference Klewicki, Fife, Wei and McMurtry2007).

To test this explanation, we plotted the Reynolds number trend of the $-{\rm \pi} /2$ phase location for an alternative wavenumber pair, $(k_x,k_z) = (1, 3)$, which was shown to more accurately identify the $y_{op}$ location for the ${Re} = 5200$ case in figure 7(b), with less than $5\,\%$ relative error. The trend for this wavenumber pair is shown by the blue squares in figure 8 and indicates that the over-saturation can be rectified by utilizing large-scale motion (LSM) modes with larger spanwise extent relative to streamwise extent. Again, this indicates that non-self-similar LSMs play an important role in accurate prediction of the phase profile and the zero-crossing location of the amplitude modulation.

Physically, we can interpret the trend of the increasing $-{\rm \pi} /2$ phase location as a reduction in the divergence rate between large and small scales (moving away from the wall, described above), or equivalently an increase in the ease of modulation between the velocity and stress modes. The fact that this phase location saturates faster for lower $k_x$ and slower for higher $k_x$ suggests that large, energetic velocity and stress modes are more easily aligned with each other as the Reynolds number increases and diffusive delays in the scale interactions become negligible. But under this interpretation, the $-{\rm \pi} /2$ phase location reflects the ease of modulation and thus we would not expect an asymptotic limit to the growth in that location with increasing Reynolds number.