2.1. Linearized Navier–Stokes equations

A statistically stationary incompressible TBL develops over a flat plate with a specified free stream velocity, $U_{\infty }(x)$. The streamwise, wall-normal and spanwise coordinates are $x$, $y$ and $z$ with $y=0$ denoting the location of the wall and $z$ assumed to be statistically homogeneous. The streamwise, wall-normal and spanwise velocities are $u$, $v$ and $w$ and the pressure divided by density is $p$. The domain is taken as a subset of the flat plate such that $x \in [x_{i},x_{o}]$, where the leading edge and virtual origin of the TBL are both outside of the domain and $y \in [0,y_{max}]$, where $y_{max}$ is in the free stream. Within the domain, reference locations, $x_{r}$, are used to define the length and velocity scales that non-dimensionalize the equations. Variables without superscripts are non-dimensionalized with the reference free stream velocity, $U_{r} = U_{\infty }(x_{r})$, and the reference length scale, $\delta _{99,r} = \delta _{99}(x_{r})$, where $\delta _{99}$ denotes the boundary layer thickness. Variables with $+$ superscripts are non-dimensionalized with the reference friction velocity $u_{\tau,r} = u_{\tau }(x_{r})$ and reference viscous length scale $\ell _{\nu } = \nu /u_{\tau,r}$, where $u_{\tau } = \sqrt {\tau _{W}/\rho }$. These reference length and velocity scales are then used to define the local outer Reynolds number, $Re = U_{r}\delta _{99,r}/\nu$ and local friction Reynolds number, $Re_{\tau } = u_{\tau,r}\delta _{99,r}/\nu$.

The instantaneous flow field, $\breve {\boldsymbol {q}}(\boldsymbol {x},t) = [\breve {\boldsymbol {u}},\breve {p}]^{T}$, can be written as the sum of a known mean flow field, $\bar {\boldsymbol {Q}}(x,y) = [\bar {\boldsymbol {U}}(x,y),\bar {P}(x,y)]^{T}$, and a fluctuation, $\boldsymbol {q}(\boldsymbol {x},t) = [\boldsymbol {u}(\boldsymbol {x},t),p(\boldsymbol {x},t)]^{T}$. For the boundary layer, $\bar {\boldsymbol {U}} = \bar {U}(x,y)\boldsymbol {e}_{x} + \bar {V}(x,y)\boldsymbol {e}_{y}$. Due to the spanwise homogeneity and statistical stationarity of the flow, the NSE may be Fourier transformed in $z$ and $t$. The Fourier transformed variables are denoted by $\hat {\cdot }$, such that

(2.1)\begin{equation} \hat{g}(x,y;k_{z},\omega) = \int_{t}\int_{z} g(x,y,z,t) \, {\rm e}^{\mathrm{i}(\omega t - k_{z} z)} \, {\rm d} z \, {\rm d} t, \end{equation}

where $k_{z} = 2{\rm \pi} /\lambda _{z}$ and $\omega = 2{\rm \pi} /\lambda _{t}$ are the spanwise wavenumber and temporal frequency. $\lambda _{z}$ and $\lambda _{t}$ then denote the spanwise and temporal wavelengths. The equations governing the Fourier modes with $(\omega,k_{z}) \ne (0,0)$ are then

(2.2)\begin{equation} \left[\begin{array}{cc} -\mathrm{i}\omega + \bar{\boldsymbol{U}} \boldsymbol{\cdot} \boldsymbol{\nabla} + \boldsymbol{\nabla}\bar{\boldsymbol{U}} -Re^{{-}1}\nabla^{2} + \epsilon_{s} & \boldsymbol{\nabla}\\ \boldsymbol{\nabla} \boldsymbol{\cdot} & 0 \end{array}\right] \hat{\boldsymbol{q}} = \left[\begin{array}{c} -\boldsymbol{\nabla} \boldsymbol{\cdot} \widehat{\boldsymbol{u}\boldsymbol{u}}\\0\end{array}\right],\end{equation}

where $\boldsymbol {\nabla } = [\partial _{x},\partial _{y},\mathrm {i} k_{z}]$, $\nabla ^{2} = \partial _{xx} + \partial _{yy} - k_{z}^{2}$ and $\epsilon _{s}$ is the sponge used for artificial boundary conditions. Following McKeon & Sharma (Reference McKeon and Sharma2010), the nonlinear terms in (2.2) are treated as a nonlinear forcing with spanwise wavenumber $k_{z}$ and temporal frequency $\omega$, $\hat {\boldsymbol {f}} = [\hat {\boldsymbol {f}}_{\boldsymbol {u}},0]^{T}$. To focus on what insight can be obtained from the linear resolvent operator, the forcing is assumed to be uncorrelated from the turbulent fluctuations. While eddy viscosity can be used to incorporate turbulent interactions into the LNSE (Del Alamo & Jimenez Reference Del Alamo and Jimenez2006; Cossu et al. Reference Cossu, Pujals and Depardon2009; Hwang & Cossu Reference Hwang and Cossu2010), their specific forms are data-driven or ad hoc and can obfuscate the true linear amplification in the NSE and the interpretability of the nonlinearities (Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021). In this approach, $\hat {\boldsymbol {f}}$ can be interpreted as an externally applied body force.

No-slip and no-penetration boundary conditions are applied such that $\hat {\boldsymbol {u}}(x,y=0) = \boldsymbol {0}$ is enforced. In the free stream, Neumann boundary conditions are applied such that $\partial _{y}\hat {\boldsymbol {u}}(x,y=y_{max}) = \boldsymbol {0}$. Since the streamwise domain is a subset of the flat plate with $x \in [x_{i},x_{o}]$, artificial boundary conditions are applied in the streamwise direction such that $\hat {\boldsymbol {q}}$ remains compact within the domain. This is enforced through the sponge, $\epsilon _{s}(x)$, in (2.2), which is $0$ for $90\,\%$ of the domain and ramps up quadratically from $0$ to $\epsilon _{0}$ at the endpoints of the domain. The explicit form of $\epsilon _s(x)$ is

(2.3) \begin{equation} \epsilon_s(x)= \left\{\begin{array}{@{}ll} \epsilon_{0}\dfrac{( | x - x_{c} - L/2 | -0.45L_{x} )^{2}}{0.0025L^{2}} & \text{if } | x - x_{c} - L/2 |\geq 0.45L,\\ 0 & \text{otherwise}, \end{array}\right. \end{equation}

where $L = x_{o} - x_{i}$ and $x_{c}$ is the centre of the domain. $\epsilon _{s}$ causes $\hat {\boldsymbol {q}}$ to decay to $\boldsymbol {0}$ in the vicinity of the boundaries. Similar strategies have been applied for linear analyses of incompressible flows in Ran et al. (Reference Ran, Zare, Hack and Jovanović2019) and Abreu et al. (Reference Abreu, Tanarro, Cavalieri, Schlatter, Vinuesa, Hanifi and Henningson2021). The inlet and outlet of the domains are then treated with Dirichlet boundary conditions such that $\hat {\boldsymbol {q}}(x = x_{i},y) = \hat {\boldsymbol {q}}(x = x_{o},y) = \boldsymbol {0}$ to enforce the compact nature of the resolvent modes. The effect of the sponge is qualitatively similar to fringe zones that are commonly used in TBL simulations to dampen the outflow and allow for the treatment of periodic boundary conditions (Chevalier et al. Reference Chevalier, Schlatter, Lundbladh and Henningson2007). In the calculations presented herein, $\epsilon _{s} = 30$. Varying this value was found to have a negligible effect on the linear amplification, provided that $\epsilon _{s}$ was large enough to dampen $\hat {\boldsymbol {q}}$ at the boundaries.

2.2. Resolvent analysis

Equation (2.2) can be written compactly as

(2.4)\begin{equation} (-\mathrm{i}\omega\mathcal{B} + \mathcal{L})\hat{\boldsymbol{q}} = \mathcal{B}\hat{\boldsymbol{f}}, \end{equation}

where $\mathcal {B}$ projects away the $p$ component as $\mathcal {B}\hat {\boldsymbol {q}} = [\hat {\boldsymbol {u}},0]^{T}$. The kinetic energy of the disturbances, integrated over the domain, is a norm induced by the inner product,

(2.5)\begin{equation} \langle\boldsymbol{a},\boldsymbol{b} \rangle = \frac{1}{L}\int_{x_{i}}^{x_{o}}\int_{0}^{y_{max}} \boldsymbol{a}^{*} \mathcal{B}\boldsymbol{b}\, {\rm d} y\, {\rm d}\kern0.7pt x. \end{equation}

By using this inner product, the adjoint of (2.2) can also be written compactly as

(2.6)\begin{equation} (\mathrm{i}\omega\mathcal{B} + \mathcal{L}^{{{\dagger}}})\tilde{\boldsymbol{q}} = \mathcal{B}\tilde{\boldsymbol{f}}, \end{equation}

where the daggers denote adjoint operators and tildes denote adjoint variables. An explicit form of the adjoint operator can be found in Gomez (Reference Gomez2024). It is assumed that $\omega$ is not an eigenvalue of $\mathcal {L}$ so that the resolvent, $\mathcal {H} = (-\mathrm {i}\omega \mathcal {B} + \mathcal {L})^{-1}\mathcal {B}$, and adjoint resolvent, $\mathcal {H}^{{\dagger} } = (\mathrm {i}\omega \mathcal {B} + \mathcal {L}^{{\dagger} })^{-1}\mathcal {B}$, exist.

In resolvent analysis, one is interested in finding the forcing inputs of $\mathcal {H}$ that lead to the largest amplification. Here, the forcing inputs and response outputs of $\mathcal {H}$ are measured with the induced norms $\| \cdot \|_{f}$ and $\| \cdot \|_{r}$ from the inner products $\langle \boldsymbol {a},\boldsymbol {b} \rangle _{f} = \langle \boldsymbol {a},\mathcal {W}_{f}\boldsymbol {b} \rangle$ and $\langle \boldsymbol {a},\boldsymbol {b} \rangle _{r} = \langle \boldsymbol {a},\mathcal {W}_{r}\boldsymbol {b} \rangle$, respectively. Here, $\mathcal {W}_{r}$ and $\mathcal {W}_{f}$ are positive definite operators that dictate the components and spatial domains of $\hat {\boldsymbol {q}}$ and $\hat {\boldsymbol {f}}$ that are used to measure the response and the forcing. The exact form of $\mathcal {W}_{r}$ and $\mathcal {W}_{f}$ will be specified for the studies considered. An example of the regions included in the norms $\| \cdot \|_{f}$ and $\| \cdot \|_{r}$ are depicted in figure 1.

Figure 1. Schematic depicting in white the integration domains for $\| \cdot \|_{f}$ (a) and $\| \cdot \|_{r}$ (b). The striped regions at the edges of the domain depict the sponge layers and the black solid curve depicts a representative boundary layer edge.

The linear amplification, $\sigma$, of $\mathcal {H}$ is here defined as the Rayleigh quotient,

(2.7) \begin{equation} \sigma^{2} =\frac{\| \mathcal{H}\hat{\boldsymbol{f}} \|_{r}}{\|\,\hat{\boldsymbol{f}} \|_{f}} =\frac{\langle\mathcal{H}\hat{\boldsymbol{f}},\mathcal{W}_{r}\mathcal{H} \hat{\boldsymbol{f}} \rangle}{\langle\,\hat{\boldsymbol{f}},\mathcal{W}_{f}\hat{\boldsymbol{f}} \rangle} =\frac{\langle\,\hat{\boldsymbol{f}},\mathcal{H}^{{{\dagger}}}\mathcal{W}_{r}\mathcal{H} \hat{\boldsymbol{f}} \rangle}{\langle\,\hat{\boldsymbol{f}},\mathcal{W}_{f}\hat{\boldsymbol{f}} \rangle}.\end{equation}

To find the largest amplification, one seeks the forcing, $\boldsymbol {\phi }$, that produces the largest $\sigma$ such that

(2.8)\begin{equation} \boldsymbol{\phi} = \operatorname*{argmax}_{\|\,\hat{\boldsymbol{f}} \|_{f}{=1}} \sigma. \end{equation}

It can be shown that (2.8) can be solved via the eigenvalue problem,

(2.9)\begin{equation} \mathcal{H}\mathcal{W}_{r}\mathcal{H}^{{{\dagger}}}\boldsymbol{\phi}_{j} = \sigma_{j}^{2}\mathcal{W}_{f}\boldsymbol{\phi}_{j}.\end{equation}

Since $\mathcal {H}\mathcal {W}_{r}\mathcal {H}^{{\dagger} }$ and $\mathcal {W}_{f}$ are positive definite operators, the eigenvalues, $\sigma _{j}^{2}$, are positive and the eigenvectors, $\boldsymbol {\phi }_{j}$, are orthonormal with respect to the inner product $\langle {\cdot },{\cdot } \rangle _{f}$ such that $\langle \boldsymbol {\phi }_{j},\boldsymbol {\phi }_{k} \rangle _{f} = \delta _{jk}$, where $\delta _{jk}$ is the Kronecker delta. The $\sigma _{j}$ are ordered such that $\sigma _{1} \ge \sigma _{2} \ge \dots \ge 0$. $\boldsymbol {\phi }_{j}$ are here identified as the forcing modes with corresponding linear gain, $\sigma _{j}$. The response modes, $\boldsymbol {\psi }_{j}$, are defined as $\boldsymbol {\psi }_{j} = \sigma _{j}^{-1}\mathcal {H}\boldsymbol {\phi }_{j}$. It can be shown that the response modes are orthonormal with respect to the inner product $\langle {\cdot },{\cdot } \rangle _{r}$ such that $\langle \boldsymbol {\psi }_{j},\boldsymbol {\psi }_{k} \rangle _{r} = \delta _{jk}$.

In this study, two sets of norms are considered. The first set is tailored to study the effect of local $Re_{\tau }$ and $\beta$ on the linear amplification in § 3. This requires that the response modes are supported at $x_{r_{1}}$ where $\beta (x_{r_{1}})$ and $Re_{\tau }(x_{r_{1}})$ are defined. If there is no spatial masking applied to the inner product such that $\langle {\cdot },{\cdot } \rangle _{r} = \langle {\cdot },{\cdot } \rangle _{f} = \langle {\cdot },{\cdot } \rangle$, then the response and forcing modes can be supported anywhere along the streamwise domain. In order to capture the modes at $x_{r_{1}}$, spatial masking is applied to the response modes so that

(2.10)\begin{equation} \langle\boldsymbol{a},\boldsymbol{b} \rangle_{r} = \frac{1}{L^+}\int_{x_{r_{1}}^+{-}L^+_{d}/2}^{x_{r_{1}}^+{+}L^+_{d}/2}\int_{0}^{y_{max}}\boldsymbol{a}^{*}\mathcal{B}\boldsymbol{b} \, {\rm d} y \, {\rm d}\kern0.7pt x^+,\end{equation}

where $L_{d}^+ = 2150$, while $\langle {\cdot },{\cdot } \rangle _{f} = \langle {\cdot },{\cdot } \rangle$. The forcing modes can extend across the entire domain allowing forcing from upstream the spatial mask to influence the linear response within the spatial mask. This helps capture non-local amplification mechanisms like those from the convective non-normality. Here $L_{d}$ is held constant in viscous units since the small-scale modes are expected to scale in inner units. For a mode to be optimal under the new inner products via the definition in (2.7), the optimal response modes must be supported in the streamwise region within $L_{d}^+/2$ of $x_{r_{1}}$. The effect of $L_{d}^+$ has been shown in the appendix of Gomez (Reference Gomez2024), where it is shown that increasing $L_{d}^+$ reduces the response of the near-wall small-scale modes at $x_{r_{1}}$. The second set of inner products is used to study the history effects in § 4. In order to only consider the effect of the pressure gradient upstream of $x_{r_{2}}$, the inner products $\langle {\cdot },{\cdot } \rangle _{r}$ and $\langle {\cdot },{\cdot } \rangle _{f}$ are masked such that

(2.11)\begin{equation} \langle\boldsymbol{a},\boldsymbol{b} \rangle_{r} = \langle\boldsymbol{a},\boldsymbol{b} \rangle_{f} = \frac{1}{L_{x}}\int_{x_{i}}^{x_{r_{2}}}\int_{0}^{y_{max}}\boldsymbol{a}^{*}\mathcal{B}\boldsymbol{b} \,{{\rm d} y} \,{{\rm d}\kern0.7pt x},\end{equation}

where $x_{i}$ is 12 $\delta _{99}(x_{r_{2}})$ upstream of $x_{r_{2}}$. This makes the resolvent modes optimal in the sense of (2.8) if they are supported upstream of $x_{r_{2}}$ and also prevents the amplification mechanisms present downstream of $x_{r_{2}}$ from contributing to the linear amplification.

It is important to note that $\boldsymbol {\psi }_{j}$, $\boldsymbol {\phi }_{j}$ and $\sigma _{j}$ are dependent on the $(\lambda _{z},\lambda _{t})$ pair, $\bar {\boldsymbol {U}}$, and choice of norms. Here $\boldsymbol {\psi }_{j}$ can be used to reconstruct the flow field as $\hat {\boldsymbol {q}} = \sum _{j}\xi _{j}\sigma _{j}\boldsymbol {\psi }_{j}$, where $\xi _{j}$ are coefficients that are found through fitting to data or appealing to the full NSE. The studies herein will primarily consider the optimal response, $\sigma _{1}\boldsymbol {\psi }_{1}$, as a rank-one approximation for the Fourier modes without any nonlinear closure. By taking an integral over $\omega$, one can also define the premultiplied amplification,

(2.12)\begin{equation} \hat{E}(y^+;k_{z}^+) = k_{z}^+\int_{\omega_{min}^+}^{\omega_{max}^+} | \sigma_{1}^+(k_{z}^+,\omega^+)\psi_{u,1}^+(x_{r_{1}}^+, y^+;k_{z}^+,\omega^+) |^{2} \, {\rm d} \omega^+.\end{equation}

Here $\hat {E}$ can be thought of as the premultiplied streamwise kinetic energy spectrum for a velocity signal under the rank-one representation, $\hat {\boldsymbol {u}} = \sigma _{1}\boldsymbol {\psi }_{\boldsymbol {u},1}$. Here, $\hat {E}(y^+;k_{z}^+)$ is a measure of the linear amplification across all scales.

2.3. Numerical discretization

The domain is discretized with $N_{x}$ equispaced points in the streamwise direction and $N_{y}$ points in the wall-normal direction. The wall-normal grid points are stretched such that half the points are below $y_{min}$. Each wall-normal grid point, $y_{k}$, is mapped from an equispaced grid point, $y'_{k}$ using $y_{k} = a y_{k}'/(b-y_{k}')$ where $a = y_{max}y_{min}/(y_{max}-2y_{min})$ and $b = 1 + a/y_{max}$. This is similar to the grids used in the supersonic linear analyses of Kamal et al. (Reference Kamal, Rigas, Lakebrink and Colonius2020) and Malik (Reference Malik1990), except here, the grid stretching is used to increase the resolution in the near-wall region of the TBLs studied. Both streamwise and wall-normal directions use a fourth-order summation by parts scheme for the differentiation operators (Mattsson & Nordström Reference Mattsson and Nordström2004). Once discretized, the operators $\mathcal {B}$, $\mathcal {L}$, $\mathcal {L}^{{\dagger} }$, $\mathcal {H}$ and $\mathcal {H}^{{\dagger} }$ become the $4N_{x}N_{y} \times 4N_{x}N_{y}$ matrices ${\boldsymbol{\mathsf{B}}}$, ${\boldsymbol{\mathsf{L}}}$, ${\boldsymbol{\mathsf{L}}}^{{\dagger} }$, ${\boldsymbol{\mathsf{H}}}$ and ${\boldsymbol{\mathsf{H}}}^{{\dagger} }$.

In order to perform the integration in (2.5), the streamwise and wall-normal integration is performed with a trapezoidal scheme. The positive definite matrix, ${\boldsymbol{\mathsf{W}}}_{T}$, performs the integration in (2.5) such that

(2.13)\begin{equation} \langle\boldsymbol{a},\boldsymbol{b} \rangle \approx \boldsymbol{a}^{*}{\boldsymbol{\mathsf{W}}}_{T}\boldsymbol{b}. \end{equation}

In a similar fashion, the operators $\mathcal {W}_{r}$ and $\mathcal {W}_{f}$ are discretized with the diagonal matrices ${\boldsymbol{\mathsf{M}}}_{r}$ and ${\boldsymbol{\mathsf{M}}}_{f}$. The diagonal entries of ${\boldsymbol{\mathsf{M}}}_{r}$ and ${\boldsymbol{\mathsf{M}}}_{f}$ are $0$ for the entries omitted and $1$ for the entries included in the inner products. As a result, the inner products for the discretized space are defined as

(2.14)\begin{equation} \langle\boldsymbol{a},\boldsymbol{b} \rangle_{r} = \boldsymbol{a}^{*}{\boldsymbol{\mathsf{M}}}_{r}{\boldsymbol{\mathsf{W}}}_{T}\boldsymbol{b} = \boldsymbol{a}^{*}{\boldsymbol{\mathsf{W}}}_{r}\boldsymbol{b}, \end{equation}

and

(2.15)\begin{equation} \langle\boldsymbol{a},\boldsymbol{b} \rangle_{f} = \boldsymbol{a}^{*}{\boldsymbol{\mathsf{M}}}_{f}{\boldsymbol{\mathsf{W}}}_{T}\boldsymbol{b} = \boldsymbol{a}^{*}{\boldsymbol{\mathsf{W}}}_{f}\boldsymbol{b}. \end{equation}

In the discretized space, the eigenvalue problem in (2.9) becomes

(2.16)\begin{equation} {\boldsymbol{\mathsf{H}}}{\boldsymbol{\mathsf{W}}}_{r}{\boldsymbol{\mathsf{H}}}^{{{\dagger}}}\boldsymbol{\phi}_{j} = \sigma_{j}^{2}{\boldsymbol{\mathsf{W}}}_{f}\boldsymbol{\phi}_{j},\end{equation}

while $\boldsymbol {\psi }_{j} = \sigma _{j}^{-1}{\boldsymbol{\mathsf{H}}}{\boldsymbol {\phi }_{j}}$. Equation (2.16) is solved using the Arnoldi algorithm (Saad Reference Saad2011). The matrices ${\boldsymbol{\mathsf{B}}}$, ${\boldsymbol{\mathsf{L}}}$, ${\boldsymbol{\mathsf{L}}}^{{\dagger} }$, ${\boldsymbol{\mathsf{W}}}_{r}$ and ${\boldsymbol{\mathsf{W}}}_{f}$ are sparse, while the resolvent matrices, ${\boldsymbol{\mathsf{H}}}$ and ${\boldsymbol{\mathsf{H}}}^{{\dagger} }$, are dense matrices that are costly to compute. The Arnoldi approach requires repeated calculations of $\boldsymbol {v} = {\boldsymbol{\mathsf{H}}}\boldsymbol {b}$, or equivalently solving $(-\mathrm {i}\omega {\boldsymbol{\mathsf{B}}} + {\boldsymbol{\mathsf{L}}})\boldsymbol {v} = {\boldsymbol{\mathsf{B}}}\boldsymbol {b}$, and similarly for the adjoint. Following the approach described in Sipp & Marquet (Reference Sipp and Marquet2013) and Schmidt et al. (Reference Schmidt, Towne, Rigas, Colonius and Brès2018), LU decompositions of $(-\mathrm {i}\omega {\boldsymbol{\mathsf{B}}} + {\boldsymbol{\mathsf{L}}})$ and $(\mathrm {i}\omega {\boldsymbol{\mathsf{B}}} + {\boldsymbol{\mathsf{L}}}^{{\dagger} })$ are computed once. The LU factors can then efficiently solve $(-\mathrm {i}\omega {\boldsymbol{\mathsf{B}}} + {\boldsymbol{\mathsf{L}}})\boldsymbol {v} = {\boldsymbol{\mathsf{B}}}\boldsymbol {b}$ through Gaussian elimination. Thus, the overall cost of the Arnoldi algorithm is the cost of the LU decompositions (Jeun et al. Reference Jeun, Nichols and Jovanović2016). The LU decomposition and Gaussian elimination is computed using the Intel oneAPI Math Kernel Library PARDISO (Schenk & Gärtner Reference Schenk and Gärtner2004), which is optimized for the sparse form of the matrices and uses parallelized subroutines. See the appendix of Gomez (Reference Gomez2024) for an algorithm explaining this strategy in more detail.

2.4. Description of datasets

This study focuses on mild to moderate APG TBL where $\beta \in [0,5]$. Since the large-scale modes require large streamwise domains to both be resolved and achieve significant amplification, only studies with large, well-resolved streamwise domains, relative to $\delta _{99}(x_{r})$, are considered. The mean flow fields of interest are the large eddy simulation flat plate TBL data of Bobke et al. (Reference Bobke, Vinuesa, Örlü and Schlatter2017), Pozuelo et al. (Reference Pozuelo, Li, Schlatter and Vinuesa2022) and Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014) which each have domains upstream of $x_{r_{2}}$ larger than $20\delta _{99}(x_{r_{2}})$. The datasets all use a constant free stream velocity, $U_{\infty }(x)$, followed by an APG imposed by the boundary condition $U_{\infty }(x) = (x-x_{0})^{m}$, where $x_{0}$ and $m$ are specified to create varying degrees of APG strengths. Under this regime, the APG TBLs are in the ‘near-equilibrium’ regime (Mellor & Gibson Reference Mellor and Gibson1966; Townsend Reference Townsend1976; Bobke et al. Reference Bobke, Vinuesa, Örlü and Schlatter2017).

A brief description of the parameters of interest in the datasets are described in table 1. In this study, the Clauser parameter is defined as

(2.17)\begin{equation} \beta(x) = \frac{\delta^{*}(x)}{u_{\tau}^{2}(x)}\frac{{\rm d} P}{{\rm d}\kern0.7pt x}(x,y=\delta_{99}(x)),\end{equation}

as in Bobke et al. (Reference Bobke, Vinuesa, Örlü and Schlatter2017) and Pozuelo et al. (Reference Pozuelo, Li, Schlatter and Vinuesa2022). Pozuelo et al. (Reference Pozuelo, Li, Schlatter and Vinuesa2022) notes that the free stream boundary conditions in the simulations cause $\bar {U}$ to decrease in the wall-normal location from its peak at the boundary layer edge such that ${\rm d}P/{{\rm d}\kern0.7pt x}$ at $y = \delta _{99}$ is different from the imposed pressure gradient (PG) through the free stream boundary condition. Due to this, the parameter

(2.18)\begin{equation} \beta_{\infty}(x) = \frac{\delta^{*}(x)}{u_{\tau}^{2}(x)}\frac{{\rm d}P}{{\rm d}\kern0.7pt x}(x,y_{max}) = \frac{\delta^{*}(x)}{u_{\tau}^{2}(x)}\frac{{\rm d}P_{\infty}}{{\rm d}\kern0.7pt x}(x),\end{equation}

will also be used to account for the effect of the free stream boundary condition. Here, $\beta _{\infty }$ is used to parameterize the effect of history by accounting for the streamwise variation of the imposed APG boundary condition. In this study, the history effect is quantified by

(2.19) \begin{align} \langle \beta_{\infty} \rangle &= \frac{1}{x_{r_{2}} - x_{i,2}}\int_{x_{i,2}}^{x_{r_{2}}} \beta_{\infty}(x)\, {{\rm d}\kern0.7pt x}\nonumber\\ &= \frac{1}{x_{r_{2}} - x_{i,2}}\int_{x_{i,2}}^{x_{r_{2}}} \int_{0}^{\infty} \frac{1}{u_{\tau}{^2}(x)}(1-\bar{U}(x,y))\frac{{\rm d}P_{\infty}}{{\rm d}\kern0.7pt x}(x) \, {\rm d} y \, {\rm d}\kern0.7pt x,\end{align}

where $x_{i,2}$ is an upstream location. Here, $x_{i,2} = x_{r_{2}} - 12\delta _{99}(x_{r_{2}})$ corresponds to the inlet location of the domains used to study the history effect in § 4. Here $\langle \beta _{\infty } \rangle$ is similar to the streamwise average that Vinuesa et al. (Reference Vinuesa, Örlü, Sanmiguel Vila, Ianiro, Discetti and Schlatter2017) considered, though they averaged $\beta$ over the momentum thickness-based Reynolds number. The domains used in Vinuesa et al. (Reference Vinuesa, Örlü, Sanmiguel Vila, Ianiro, Discetti and Schlatter2017) for the averaging were $O(10\delta _{99})$ large, which is consistent with the averaging domain used here. Furthermore, (2.19) offers a measure of the APG history in a manner analogous to the streamwise averaging induced by the inner product in (2.11).

Table 1. Parameters describing the datasets, along with quantities of interest. The parameters $m$ and $x_{0}$ set the PG in each simulation. Note that the $Re_{\tau }$ range is listed for the entire dataset, although only subsets of the domains are considered in this study. The streamwise locations $x_{r_{1}}$ and $x_{r_{2}}$ are chosen such that $Re_{\tau }(x_{r_{1}}) = 550$ and $Re_{\tau }(x_{r_{2}}) = 777$. Here $\langle \beta _{\infty } \rangle$ is computed following (2.19). The ZPG dataset, S0, is from Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014) and the APG dataset, b14, is from Pozuelo et al. (Reference Pozuelo, Li, Schlatter and Vinuesa2022) while the rest are from Bobke et al. (Reference Bobke, Vinuesa, Örlü and Schlatter2017).

Figure 2(a) shows the variation of $\beta$ with $Re_{\tau }$ for the datasets used, along with the values at the reference locations used in this study. Figure 2(b–d) show $\bar {U}^+$, $\overline {uu}^+$ and $\bar {V}^+$ at fixed $x_{r_2}$. Consistent with the observations from previous studies described above, these mean flow fields demonstrate monotonic growth in the wake velocity and decreasing $\bar {U}^+$ in the log-layer region with increasing $\beta$. Monotonic growth in the secondary peak of $\overline {uu}^+$ is observed with $\langle \beta _{\infty } \rangle$, with the near-wall peak of $\overline {uu}^+$ also affected. The non-monotonic behaviour in $\bar {V}^+$ with $\langle \beta \rangle _{\infty }$ highlight the non-equilibrium nature of the datasets, with different levels of deceleration in the flow at $x_{r_{2}}$

Figure 2. (a) Here $\beta$ plotted against $Re_{\tau }$ for the APG datasets described in table 1. The solid lines indicate the domain used to study the history effect. The values evaluated at $x_{r_{2}}$ and $x_{r_{1}}$ are plotted with circles and squares, respectively. The thicker lines denote the subdomains used in § 4 to study the history effect. Plots of $\bar {U}^+$ (b), $\overline {uu}^+$ (c) and $\bar {V}^+$ (d) against $y^+$ at $x_{r_{2}}$. The viscous units are also computed at $x_{r_{2}}$.

3.1. Effect of masking on the response and forcing modes

Here, the goal is to study the linear response of large- and small-scale modes subject to different local parameters. This requires that the response modes are supported at $x_{r_{1}}$ where $\beta (x_{r_{1}})$ and $Re_{\tau }(x_{r_{1}})$ are defined. Without a spatial mask, the response and forcing modes can be supported anywhere along the streamwise domain. It was shown that the small-scale modes, characterized by $O({\lambda _{z}^+}) \le O({100})$, were amplified via the componentwise non-normality from the mean shear near the wall while the large-scale modes, characterized by $O({\lambda _{z}}) \ge O({1})$, extended into the outer region of the flow and were amplified via the convective non-normality (Gomez Reference Gomez2024). As a result, the unmasked small-scales are expected to be concentrated in the upstream region where they experience the largest mean shear within the domain. The unmasked large-scale modes can extend across the entire domain, with the largest forcing amplitude concentrated upstream and the largest response amplitude concentrated downstream (Sipp & Marquet Reference Sipp and Marquet2013; Symon et al. Reference Symon, Rosenberg, Dawson and McKeon2018).

As an example for how the leading resolvent modes behave without masking, $\boldsymbol {\psi }_{1}$ and $\boldsymbol {\phi }_{1}$ are computed for a representative small-scale mode with $(\lambda _{z}^+,\lambda _{t}^+) = (50,50)$ and a representative large-scale mode with $(\lambda _{z},\lambda _{t}) = (1,1.4)$ using case 3. In figure 3(a,b), $\phi _{v,1}$ and $\sigma _{1}\psi _{u,1}$ are plotted for the representative small-scale and in figure 3(c,d) for the representative large-scale. As alluded to earlier, without masking, the small-scale modes are concentrated upstream of $x_{r_{1}}$ while the large-scale mode is supported at $x_{r_{1}}$. For these unmasked modes, the amplitude of $\psi _{1}(y,x_{r_{1}})$ is negligible for the small-scale mode and serves as a poor measure of the linear response at $x_{r_{1}}$.

Figure 3. Real parts of $\phi _{v,1}$ (a,c,e,g) and $\sigma _{1}\psi _{u,1}$ (b,d,f,h) for the representative small-scale (a,b,e,f) and large-scale (c,d,g,h) modes using case 3. The modes in (a–d) are computed without masking in the inner product while the modes in (e–h) are computed using masking from (2.10) in the inner product for the response modes. The edges of the spatial mask are plotted with the black solid vertical lines. The dotted line is the reference location, $x_{r_{1}}$. The curved solid line denotes the local boundary layer thickness.

In figure 3(e,f), the modes from figure 3(a,b) are recomputed using the spatial masking described in (2.10). Both the representative small-scale $\phi _{v,1}$ and $\psi _{u,1}$ are supported at $x_{r_{1}}$ with the spatial mask. Due to the larger $Re_{\tau }$ at $x_{r_{1}}$ compared with the inlet location and fixed $\lambda _{z}$ and $\lambda _{t}$ across the domain in the calculation, the modes plotted in figure 3(e,f) have significantly smaller wall-normal extents and streamwise wavelengths than those in figure 3(a,b). This is because the small-scale modes are amplified by the large mean shear concentrated in the viscous subregion, whose wall-normal extent scales with $\ell _{\nu }$. As $Re_{\tau }$ increases downstream, $\ell _{\nu }$ decreases and constrains the length scales of the modes. Now the representative large-scale modes from figure 3(c,d) are recomputed using the same spatial masking and plotted in figure 3(g,h). Since $\langle {\cdot },{\cdot } \rangle _{f}$ is left unchanged, the modes can still be amplified via the convective non-normality. As a result, the large-scale forcing can still be supported upstream of $x_{r_{1}}$ while its response is supported downstream of $x_{r_{1}}$. However, because of the spatial masking, the optimal forcing is primarily concentrated upstream of $x_{r_{1}}$ since any forcing downstream of $x_{r_{1}}+L_{d}/2$ would lead to a response outside of the spatial mask without contributing to $\sigma _{1}$. This also highlights the non-local nature in the amplification of large-scale structures and their susceptibility to upstream history effects that is explored in more detail in § 4.

For the rest of this section, the spatial masking described in (2.10) is applied to identify optimal small and large-scale responses at $x_{r_{1}}$.

3.2. Effect of $\beta$ on linear amplification at low $Re_{\tau }$

To investigate the effect of local $\beta (x_{r_{1}})$, cases 1–4 are examined since they have a similar $Re_{\tau }(x_{r_{1}})$ and different values of $\beta (x_{r_{1}})$. Furthermore, the slow variation in $\beta (x)$ near $x_{r_{1}}$ minimizes the influence of ${{\rm d}\beta }/{{\rm d}\kern0.7pt x}$ or history on the linear amplification between cases $1$–$4$. Here, the effect of $\beta (x_{r_{1}})$ on the resolvent amplification is compared for a wide range of scales.

First, the streamwise component of the linear response of the ZPG TBL of case 1 and the APG TBL of case 3 are compared for a large-scale mode ($\lambda _{t}^+ = 200, \lambda _{z} = 1$) and small-scale mode ($\lambda _{t}^+ = 50, \lambda _{z}^+ = 50$) in figure 4. Plotted in figure 4(a–c), the effect of the APG on the large-scale response at $x_{r_{1}}$ is an increase in the magnitude and peak amplitude location farther from the wall than the ZPG counterparts. Since these large-scale modes have support, i.e. non-negligible amplitude, in the outer region of the TBL, it is expected that they are most susceptible to the APG effects. On the other hand, the near-wall small-scale modes plotted in figure 4(d–f) are only slightly affected by the change in the mean flow field. This is because for this mild APG, the near-wall is still dominated by viscous effects and the APG hardly affects the inner scaling of $\bar {U}^+$. These near-wall resolvent modes have been shown to be self-similar in the region where the mean flow field demonstrates self-similarity (Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013; Gomez Reference Gomez2024).

Figure 4. Real parts of $\sigma _{1}^+\psi ^+_{u,1}$ (a,b,d,e) for large-scale (a,b) and small-scale (d,e) modes using case 1 (a,d) and case 3 (b,e). The black lines are the same as those plotted in figure 3. Line plots of $| \sigma _{1}^+\psi ^+_{u,1}(y^+,x_{r_{1}}) |$ for the large-scale (c) and small-scale (f) modes. Black and blue denote modes computed using cases 1 and 3, respectively.

While it is illustrative to compare the effect of $\beta$ on individual modes, we now consider the effect of increasing $\beta$ over a range of scales. The leading masked response modes are computed over a range of $31$ logarithmically spaced $k_{z}^+$ between $2{\rm \pi} /9.5$ and $2{\rm \pi} /11\,000$ and $33$ logarithmically spaced $\omega ^+$ between $2{\rm \pi} /1.25$ and $2{\rm \pi} /20\,000$. This range of $k_{z}^+$ and $\omega ^+$ is wide enough to capture the small-scale and large-scale modes and encompasses the characteristic length and time scales of the expected energy-containing structures in the actual TBLs. The premultiplied amplification spectra, $\hat {E}(y^+,k_{z}^+)$, (see (2.12)) are then computed for each $\beta$.

In figure 5, $\hat {E}^+$ is computed for cases 1–4. For the ZPG in case 1, $\hat {E}^+$ identifies only a near-wall peak with a local maximum near $y^+ = 20$ and $\lambda _{z}^+ = 77$. The lack of a secondary peak in the premultiplied amplification at this low $Re_{\tau }$ is similar to what is observed in the true premultiplied kinetic energy spectra for low $Re_{\tau }$ ZPG TBL and channels, which are dominated by the near-wall cycle. For the APG TBLs in cases 2-4, a near-wall peak is also observed. For case 2, this local maximum occurs at the same location as case 1, while for cases 3 and 4, the local maxima is at $y^+ = 20$ and $\lambda _{z}^+= 60$.

Figure 5. Here $\hat {E}^+(y^+;k_{z}^+)$ is plotted against $y^+$ and $\lambda _{z}^+$ for case 1 (a), case 2 (b), case 3 (c) and case 4 (d); the value of $\beta (x_{r_{1}})$ increases monotonically from (a) to (d). The black contour lines in each subplot are spaced as $95\,\%$, $82\,\%$, $50\,\%$, $10\,\%$, $1\,\%$ and $0.1{\%}$ of the maximum value of $\hat {E}^+(y^+;\lambda _{z}^+)$.

Here $\hat {E}^+$ predicts higher wall-normal locations and smaller wavelengths for the near-wall peak than observed in the turbulent data of ZPG and mild to moderate APG TBLs. Besides the assumption of uncorrelated (uncoloured) forcing, one reason for this discrepancy is that in cases 1–4, $\delta _{99}$ varies substantially while $\lambda _{z}$ remains fixed across the domain. As an example, in case 3, the $\lambda _{z}^+ = 100$ mode has a $\lambda _{z} = 0.19\delta _{99}(x_{r_{1}}) = 1.20\delta _{99}(x_{i})$. That is, near the inlet, the $\lambda _{z}^+ = 100$ mode is larger than the $\delta _{99}$ at the inlet. This in turn causes the $\lambda _{z}^+=100$ mode at $x_{r_{1}}$ to behave like a large-scale mode and experiences amplification via the convective non-normality, rather than the componentwise non-normality characteristic of the small-scale modes. As a result, the $\lambda _{z}^+$ modes are absent from the near-wall peak for these low $Re_{\tau }$ TBLs. The changing nature of the mode behaviour with streamwise distance is a challenge unique to the TBL relative to the canonical internal flows. It was shown in Gomez (Reference Gomez2024) that by also applying masking to the forcing such that $\langle {\cdot },{\cdot } \rangle _{f}= \langle {\cdot },{\cdot } \rangle _{r}$, a local maximum could be found in $\hat {E}^+$ with $\lambda _{z}^+ \approx 100$ for both ZPG and APG TBLs.

Despite the simplified forcing assumption, $\hat {E}^+$ also demonstrates a secondary peak for large-scale modes farther from the wall for the APG TBLs in cases 2–4. The secondary peak in $\hat {E}^+$ is amplified monotonically with increasing $\beta$. Similarly, a monotonic increase in the secondary peak of the spanwise kinetic energy spectra was observed in the DNS data of Lee (Reference Lee2017) at similar $Re_{\tau }$ to the mean data used here. For larger $Re_{\tau }$, the streamwise energy spectra in the experiment also demonstrate large-scale secondary peaks that are energized with increasing $\beta$ (Sanmiguel Vila et al. Reference Sanmiguel Vila, Vinuesa, Discetti, Ianiro, Schlatter and Örlü2020). Despite the similar trends seen with $\beta$, there is again a discrepancy in the location of the secondary peaks in $\hat {E}^+$ and those of similar low $Re_{\tau }$ spanwise energy spectra (Bobke et al. Reference Bobke, Vinuesa, Örlü and Schlatter2017; Lee Reference Lee2017). The secondary peaks in figure 5(b–d) occur at $\lambda _{z}^+ \approx 0.5Re_{\tau }$ whereas the turbulent data reports $\lambda _z^+ \approx 0.8Re_{\tau }-0.9Re_{\tau }$. The discrepancy is likely attributable to the lack of nonlinear closure, but may again be in part related to the change in $\delta _{99}$ across the domain, as a mode of the order of $\delta _{99}(x_{r_{1}})$ at $x_{r_{1}}$ would be around $6\delta _{99}(x_{i})$ for case 3. Since the large-scale modes are forced from the upstream location, the smaller upstream boundary layer can affect the behaviour of the amplification at $x_{r_{1}}$ by geometrically constraining the modes that can be amplified downstream. Nonetheless, the masked biglobal resolvent approach is able to capture the presence of both a small-scale, near-wall peak and a large-scale peak that is amplified monotonically with increasing $\beta$.

3.3. Effect of $Re_{\tau }$ on linear amplification

Linear amplification is now considered using the ZPG TBL dataset of Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014) for the low $Re_{\tau }$ case 1 and larger $Re_{\tau }$ case 6. The large $Re_{\tau }$ APG TBL dataset of Pozuelo et al. (Reference Pozuelo, Li, Schlatter and Vinuesa2022) (case 7) is also used to study the $\beta$ effects in the linear amplification at higher $Re_{\tau }$ by comparing with case 5 at the matched $Re_{\tau } = 1500$. The sweeps are computed using the same $k_{z}$ and $\omega$ range as the previous results while the wall-normal stretching has been adjusted to resolve the smaller inner region in the larger $Re_{\tau }$ datasets by shrinking $y_{min}$.

Figure 6(a,b) show the variation of $\hat {E}^+$ for cases 1 and 6. Along with the near-wall peaks, $\hat {E}^+$ also demonstrate increased amplification for the large-scale modes with increasing $Re_{\tau }$ with case 6 showing evidence of a secondary peak. This secondary peak occurs at a wall-normal location slightly closer to the wall than the secondary peaks of the APG results in figure 5, potentially pointing to a difference between log and APG-amplified wake structures. However, the $Re_{\tau }$ is too small to adequately distinguish between the two. Figure 6(d) provides direct comparison of isocontours of $\hat {E}^+$ for cases 1 and 6. These plots demonstrate that $\hat {E}^+$ is self-similar in the near-wall region due to the inner-scaled mask in $\langle {\cdot },{\cdot } \rangle _{r}$ and self-similar mean flows, plotted in figure 6(c) for reference. The near-wall self-similarity of the resolvent amplification in individual modes had been studied for channel flows in Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013) and spatially developing TBLs in Gomez (Reference Gomez2024) where the response and forcing modes are spatially constrained via small, inner-scaled, streamwise domains. Figure 6(d) demonstrates that the self-similarity in the resolvent modes persists in $\hat {E}^+$ even when the forcing modes are not spatially constrained with a spatial mask or limited domain length. As the large-scale outer-scaled resolvent modes become more amplified farther from the wall, the near-wall self-similarity of $\hat {E}^+$ breaks because the inner-scaled self-similarity of $\bar {\boldsymbol {U}}$ no longer holds.

Figure 6. Here $\hat {E}^+(y^+;k_{z}^+)$ is plotted against $y^+$ and $\lambda _{z}^+$ for case 1 (a) and case 6 (b). Here $\bar {U}^+(x_{r_{1}},y^+)$ (c) and contour lines of $\hat {E}^+(y^+;k_{z}^+)$ (d) for cases 1 and 6 plotted in dotted and solid lines, respectively. The contours in (a), (b) and (d) denote $95\,\%$, $82\,\%$, $50\,\%$, $10\,\%$, $1\,\%$ and $0.1{\%}$ of the maximum value of $\hat {E}_{uu}^+(y^+;\lambda _{z}^+)$.

Now, the effect of $\beta$ is considered for the large $Re_{\tau }$ APG TBL from case 7. These results are directly compared with those of case 5 which is at the same $Re_{\tau }$, with their mean velocity profiles compared in figure 7(a). At this larger $Re_{\tau }$ and mild $\beta$, the inner region of the mean flow is hardly affected by the APG. As a result, $\hat {E}^+$ exhibits near-wall self-similarity in figure 7(b), like case 6. For $\hat {E}^+$, the self-similarity breaks at larger $\lambda _{z}^+$, indicating that these large-scale modes are more sensitive to the changes in the mean flow field than the near-wall small-scale modes. Just as was seen in figure 5, the large-scale modes are more amplified for the APG TBL than the ZPG TBL in the higher $Re_{\tau }$ case, similar to the energization of the outer-scaled spectral peak in the APG TBL observed in Pozuelo et al. (Reference Pozuelo, Li, Schlatter and Vinuesa2022).

Figure 7. Here $\bar {U}^+(x_{r_{1}},y^+)$ (a) and contour lines of $\hat {E}^+(y^+;k_{z}^+)$ (b) for cases 5 and 7 plotted in black and orange lines, respectively. The contours in (b) denote $95\,\%$, $82\,\%$, $50\,\%$, $10\,\%$, $1\,\%$ and $0.1$ of the maximum value of $\hat {E}_{uu}^+(y^+;\lambda _{z}^+)$.

The masked biglobal resolvent analysis approach is able to capture $Re_{\tau }$ trends seen in the spanwise energy spectra. A self-similar near-wall peak and increased amplification in the large-scale structures with $Re_{\tau }$ are observed. At matched $Re_{\tau }$, the linear amplification increases with $\beta$. For sufficiently large $Re_{\tau }$, a secondary peak is also observed. Although this analysis is is restricted to a rank-one approximation, it nonetheless suggests that changes in the linear amplification from a change in the mean flow field due to $Re_{\tau }$ and $\beta$ can explain features seen in experimental and simulation data without any nonlinear closure in the form of a scale-dependent coefficient. It can thus be surmised that these are phenomena associated with the dominant behaviour of the linear resolvent whose predictions can be improved by more sophisticated modelling of the details of the nonlinear forcing.

5.1. Neglecting the non-parallel terms to explain the linear growth

By neglecting $\bar {V}$, $\bar {V}_{x}$ and $\bar {U}_{x}$, the LNSE in (2.2) can be approximated using the simplified forced biglobal OSS system,

(5.2)\begin{gather} \hat{v} = \mathcal{L}_{OS}^{{-}1}\,\hat{f}_{v} , \end{gather}

(5.3)\begin{gather}\hat{\omega}_{2} = \mathcal{L}_{SQ}^{{-}1}( \,\hat{f}_{\omega} + \mathrm{i} k_{z}\bar{U}_{y}\mathcal{L}_{OS}^{{-}1}\,\hat{f}_{v} ), \end{gather}

where $\mathcal {L}_{SQ} = \!-\!\mathrm {i}\omega + \bar {U}\partial _{x} - Re^{-1}\nabla ^{2}$, $\mathcal {L}_{OS} = (-\mathrm {i}\omega + \bar {U}\partial _{x} - ({1}/{Re})\nabla ^{2}) -2\nabla ^{-2}\partial _{y}(\bar {U}_{y}\partial _{x})$, $\hat {f}_{\omega } = \mathrm {i} k_{z}\hat {f}_{x} - \partial _{x}\hat {f}_{z}$ and $\hat {\omega }_{2} = \mathrm {i} k_{z}\hat {u} - \partial _{x}\hat {w}$. See the Appendix for more details. The resolvent gains of this system will then be denoted as $\tilde {\sigma }_{j}$. The forced biglobal OSS system has been studied by Ran et al. (Reference Ran, Zare, Hack and Jovanović2019) to identify the stochastic receptivity of a laminar boundary layer using the LNSE. In that case, terms involving $\bar {U}_{xx}$ and $\bar {V}_{x}$ were neglected, consistent with the Blasius similarity solution. Equations (5.2) and (5.3) make a stronger assumption by neglecting the non-parallel shear terms and wall-normal advection all together.

The next approximation comes from the observation that the most linearly amplified modes have small wall-normal components relative to the spanwise and streamwise components. This is plotted in figure 12 where $\| \psi _{v,1} \|_{r}$ contributes less than $5\,\%$ of $\| \boldsymbol {\psi }_{1} \|_{r}$. As a result, the contribution from (5.2) is negligible in the numerator of (2.7) and $\| \boldsymbol {\psi }_{1} \|_{r}^{2} \approx \| \boldsymbol {\psi }_{u,1} \|_{r}^{2} + \| \boldsymbol {\psi }_{w,1} \|_{r}^{2} = \| (\partial _{xx} - k_{z}^{2})\boldsymbol {\psi }_{\omega,1} \|_{r}^{2}$, where $\psi _{\omega,1} = \mathrm {i} k_{z}\psi _{u,1} - \partial _{x}\psi _{w,1}$ (see the Appendix). Because it was assumed that $k_{z}$ and the streamwise derivatives are $O(1)$, it can be further assumed that $\| \boldsymbol {\psi }_{1} \|_{r}^{2} \sim \| \boldsymbol {\psi }_{\omega,1} \|_{r}^{2}$. By using the approximations that $\psi _{v,1}$ is negligible and that only outer region forcing is optimal, the linear amplification can be approximated as

(5.4)\begin{equation} \tilde{\sigma}_{1} \approx \max_{\|\,\hat{\boldsymbol{f}} \|_{o} = 1}{\| \hat{\omega}_{2} \|_{o}}. \end{equation}

Through application of operator norm inequalities and the triangle inequality using $\| \hat {\omega }_{2} \|_{o}$ and (5.3), the approximation $\tilde {\sigma }_{1}$ is bounded above as

(5.5)\begin{equation} \tilde{\sigma}_{1} \le \| \mathcal{L}_{SQ}^{{-}1} \|_{o} ( \| \,\hat{f}_{\omega} \|_{o} + k_{z} \| \bar{U}_{y} \|_{o}\| \mathcal{L}_{OS}^{{-}1}\;\widehat{f_{v}} \|_{o}).\end{equation}

Equation (5.5) gives an upper bound to the approximation of $\tilde {\sigma }_{1}$ based on linear growth with $\bar {U}_{y}$. Similarly, if the forcing is constrained to only wall-normal forcing, it must be the case that

(5.6)\begin{equation} \max_{\|\,\hat{\boldsymbol{f}} \|_{o} = 1, \,\hat{f}_{u} = \,\hat{f}_{w} = 0}{\| \hat{\omega}_{2} \|_{o}} \le \tilde{\sigma}_{1}.\end{equation}

In other words, the largest possible gain that stems from solely wall-normal forcing is less than the optimal gain. Equation (5.6) can be written explicitly by setting $\hat {f}_{\omega } = 0$ in (5.3) to establish the lower bound as

(5.7)\begin{equation} k_{z}\| \mathcal{L}_{SQ}^{{-}1} \|_{o}\| \bar{U}_{y} \|_{o} \| \mathcal{L}_{OS}^{{-}1}\,\hat{f}_{v} \|_{o} \le \tilde{\sigma}_{1}.\end{equation}

Together, (5.5) and (5.7) suggest that $\tilde {\sigma }_{1}$ is bounded below by a quantity that scales with $\| \bar {U}_{y} \|_{o}$ and bounded above by a quantity that grows linearly with $\| \bar {U}_{y} \|_{o}$.

Figure 12. Plots of each of the terms in $\| \boldsymbol {\psi }_{1} \|_{r}^{2}$ for the representative mode from figure 8.

It must be noted that the $\bar {U}_{y}$ dependence in both (5.6) and (5.7) also depends on the amplification from $\mathcal {L}_{OS}^{-1}$ and $\mathcal {L}_{SQ}^{-1}$, measured in the outer region of the flow. To show that the amplification from $\mathcal {L}^{-1}_{SQ}$ will not produce a significant additional $\bar {\boldsymbol {U}}$ dependence on the upper and lower bounds, its resolvent gains, $\sigma _{j,SQ}$, are plotted in figure 13(a) for the representative mode from figure 8. Recall that $\sigma _{1,SQ} = \| \mathcal {L}_{SQ}^{-1} \| \approx \| \mathcal {L}_{SQ}^{-1} \|_{o}$. Here $\sigma _{j,SQ}$ has a slight $\langle \beta _{\infty } \rangle$ dependence for the optimal modes, where $\sigma _{1,SQ}$ grows by $18\,\%$ between $\langle \beta _{\infty } \rangle = 0$ to $\langle \beta _{\infty } \rangle = 3.47$. For reference, $\sigma _{1}$ grows by $187\,\%$ over the same $\langle \beta _{\infty } \rangle$ range. For the higher-order modes, the $\sigma _{j,SQ}$ variation with $\langle \beta _{\infty } \rangle$ is negligible compared with the variation in $\sigma _{5}$. Furthermore, $\mathcal {L}^{-1}_{SQ}$ is not low-rank as the first five $\sigma _{j,SQ}$ are all within the same order of magnitude. This means that it is more than likely that $(\,\hat {f}_{\omega } + \mathrm {i} k_{z}\bar {U}_{y}\mathcal {L}_{OS}^{-1}\;\hat {f}_{v})$ in (5.3) is composed of higher-order forcing modes of $\mathcal {L}^{-1}_{SQ}$ that have weak $\bar {\boldsymbol {U}}$ dependence. As a result, it can be concluded that the upper bound in (5.5) is weakly dependent on the $\bar {\boldsymbol {U}}$-dependence on $\| \mathcal {L}_{SQ}^{-1} \|_{o}$.

Figure 13. Plots of the first five resolvent gains of $\mathcal {L}_{SQ}$ (a), $\mathcal {L}_{OS}$ (b) and the full operator with componentwise masking to consider only $\hat {f}_{v}$ forcing and $\hat {u}$ and $\hat {w}$ response (c) for the representative mode from figure 8. Panels (a) and (b) are plotted against $\langle \beta _{\infty } \rangle$ while (c) is plotted against $\| \bar {U}_{y} \|_{o}$. The solid lines in (c) denote lines of best fit for each of the $\sigma _{j,v\rightarrow (u,w)}$. Note that the lighter colours denote increasing $j$.

To estimate the amplification from $\mathcal {L}^{-1}_{OS}$, componentwise masking is applied to consider the amplification from $\hat {f}_{v}$ to $\hat {v}$ so as to only consider the $\mathcal {L}_{OS}^{-1}\;\hat {f}_{v}$ term in (5.2). These resolvent gains are denoted as $\sigma _{j,OS}$ and are plotted in figure 13(b). Here $\sigma _{j,OS}$ demonstrates a slight decrease with $\langle \beta _{\infty } \rangle$ with a decrease of approximately $20\,\%$ between the $\langle \beta _{\infty } \rangle$ range. The $\sigma _{j,OS}$ also demonstrate a weak $\bar {\boldsymbol {U}}$ dependence on its leading gains. Although the amplification has contributions from $\mathcal {L}_{SQ}^{-1}$ and $\mathcal {L}_{OS}^{-1}$, it is expected that the $\bar {\boldsymbol {U}}$-dependence on the amplification only enters the $\bar {U}_{y}$ term.

To show that the $\bar {\boldsymbol {U}}$ dependence on $\mathcal {L}_{SQ}^{-1}$ and $\mathcal {L}_{OS}^{-1}$ does not affect the lower bound, figure 13(c) plots $\sigma _{1,v\rightarrow (u,w)}$, the resolvent gains when componentwise masking is applied to consider only forcing from $\hat {f}_{v}$ to a response with $\hat {u}$ and $\hat {w}$. This mimics the conditions in (5.7) which captures only amplification from the lift-up mechanism. Here, it is found that despite the $\langle \beta _{\infty } \rangle$ variation in the amplification of $\mathcal {L}_{SQ}^{-1}$ and $\mathcal {L}_{OS}^{-1}$, the growth in $\sigma _{1,v\rightarrow (u,w)}$ is linear with $\| \bar {U}_{y} \|_{o}$. Fitting a line to $\sigma _{1,v\rightarrow (u,w)}$ as

(5.8)\begin{equation} \sigma_{j,v\rightarrow(u,w)} = T_{j}\| \bar{U}_{y} \|_{o},\end{equation}

where $T_{j}$ is a constant for each $j$ in figure 13(c), shows that the approximations made in this section can explain the amplification from the lift-up mechanism. Furthermore, it demonstrates that any $\bar {\boldsymbol {U}}$ dependence on $\mathcal {L}_{OS}^{-1}$ and $\mathcal {L}_{SQ}^{-1}$ does not affect the linear growth with $\| \bar {U}_{y} \|_{o}$. The differences in the mode shapes are not accounted for in (5.8), which may explain why the lines of best fit do not cross $\sigma _{1,v\rightarrow (u,w)}$ for the ZPG TBL.

The linear fit in (5.8) holds when $\hat {f}_{\omega } = 0$ such that the assumption in (5.7) holds. For a more general case, a new linear fit is proposed based on the upper bound in (5.5) such that

(5.9)\begin{equation} \sigma_{j}(\lambda_{z},\lambda_{t}) = a_{j,0}(\lambda_{z},\lambda_{t}) + a_{j,1}(\lambda_{z},\lambda_{t})\| \bar{U}_{y} \|_{o},\end{equation}

where $a_{j,i}(\lambda _{z},\lambda _{t})$ are scale-dependent constants for each $j$. The coefficients $a_{j,0}$ and $a_{j,1}$ are analogous to $\| \mathcal {L}^{-1}_{SQ} \|_{o}\| \;\hat {f}_{\omega } \|_{o}$ and $k_{z}\| \mathcal {L}^{-1}_{SQ} \|_{o}\| \mathcal {L}^{-1}_{OS}\;\hat {f}_{v} \|_{o}$, respectively. These linear fits are shown in figure 14 for the same $\sigma _{j}(\lambda _{z},\lambda _{t})$ plotted in figure 9 demonstrating that there is indeed linear growth in $\sigma _{j}$ with $\| \bar {U}_{y} \|_{o}$ for the APG TBLs as suggested by the upper bound. Furthermore, figure 12(c) shows that the APG TBL datasets used here have mean flow fields such that

(5.10)\begin{equation} \| \bar{U}_{y} \|_{o} = d_{0} + d_{1}\langle \beta_{\infty} \rangle,\end{equation}

where $d_{0}$ and $d_{1}$ were constants. Taking both (5.9) and (5.10) recovers the linear form in $\sigma _{j}$ with $\langle \beta _{\infty } \rangle$ from (4.1). Hence, it can be concluded that the linear growth of $\sigma _{j}$ in $\| \bar {U}_{y} \|_{o}$ for the large-scale structures can be predicted from the LNSE if the large-scale modes are supported in the outer region of the flow. Equation (5.9) stems from the dominant terms of the LNSE while the linear fit in (4.1) requires that (5.10) holds, which is not guaranteed for general APG TBLs. Although the linear fit for the APG TBL modes was motivated by neglecting the non-parallel terms by considering $\tilde {\sigma }_{1}$; $\sigma _{1}$ is significantly larger than $\tilde {\sigma }_{1}$. This indicates that the non-parallel terms are non-negligible. Next, the influence of the non-parallel terms on the amplification will be examined.

Figure 14. The same as figure 9, except that the $\sigma _{j}$ are plotted against $\| \bar {U}_{y} \|_{o}$ and the lines of best fit are fitted using $\| \bar {U}_{y} \|_{o}$ for the data with $\langle \beta _{\infty } \rangle$.

5.2. Influence of the non-parallel terms on the linear amplification

Although artificially neglecting the non-parallel terms in the LNSE can still explain the linear growth with $\langle \beta _{\infty } \rangle$, figure 11(f) demonstrates that the non-parallel terms are a non-negligible source of amplification. Here, the effect of the non-parallel terms will be examined to discuss their role in the linear amplification. Equations (5.2) and (5.3) are now reconsidered, except that the $\mathcal {L}_{SQ}$ operators are augmented to incorporate the $\bar {U}_{x}$ terms as $\tilde {\mathcal {L}}_{SQ} = \mathcal {L}_{SQ} + \bar {U}_{x}$ while the $\mathcal {L}_{OS}$ operators are left unchanged because their amplification was shown to be approximately independent of $\langle \beta _{\infty } \rangle$ when using the full LNSE in figure 13(b). Incorporating $\bar {U}_{x}$ also introduces off-diagonal terms in the full OSS, though they are negligible relative to the mean shear, $k_{z}\| \bar {U}_{y} \|_{o}$. Since only $\mathcal {L}_{SQ}$ has changed, an upper and lower bound for $\sigma _{1}$ can be constructed by following the steps in § 5.1 as

(5.11)\begin{equation} k_{z}\| \tilde{\mathcal{L}}_{SQ}^{{-}1}\bar{U}_{y}\mathcal{L}_{OS}^{{-}1}\,\hat{f}_{v} \|_{o} \le \sigma_{1} \le \| \tilde{\mathcal{L}}_{SQ}^{{-}1} \|_{o} ( \| \,\hat{f}_{\omega} \|_{o} + k_{z} \| \bar{U}_{y} \|_{o}\| \mathcal{L}_{OS}^{{-}1} \;\widehat{f_{v}} \|_{o}). \end{equation}

The effect of $\bar {U}_{x}$ is to augment the amplification of $\tilde {\mathcal {L}}_{SQ}$. This can be seen by computing the resolvent amplification of $\tilde {\mathcal {L}}_{SQ}$, $\tilde {\sigma }_{1,SQ}$ and comparing with the $\sigma _{j,SQ}$ of $\mathcal {L}_{SQ}$ in figure 15(a). It can be seen that $\tilde {\sigma }_{j,SQ}>\sigma _{j,SQ}$ for $j = 1,2,3,4,5$, indicating that the inclusion of $\bar {U}_{x}$ increases the amplification by $10\,\%$. This increase in amplification then manifests itself as an increase in the amplification via the lift-up effect. To show this, $\sigma _{j,(v \rightarrow u,w)}$ is compared with $\tilde {\sigma }_{j,(v \rightarrow u,w)}$, the resolvent amplification from $\hat {f}_{v}$ to $(\hat {u},\hat {w})$ for the LNSE neglecting $\bar {U}_{x},$ $\bar {V}_{x}$ and $\bar {V}$, in figure 15(b), demonstrating that $\sigma _{j,(v \rightarrow u,w)}>\tilde {\sigma }_{j,(v \rightarrow u,w)}$. From the lines of best fit, the rate of growth with $\| \bar {U}_{y} \|_{o}$ is larger for the full LNSE ($\sigma _{j,(v \rightarrow u,w)}$) than by neglecting the non-parallel terms ($\tilde {\sigma }_{j,(v \rightarrow u,w)}$). From the lower bound in (5.11), the increased rate of growth stems from the increase in amplification of $\tilde {\mathcal {L}}_{SQ}^{-1}$ since it was argued that $\mathcal {L}_{OS}^{-1}$ was hardly affected by the change in base flow. This increased rate of growth with $\| \bar {U}_{y} \|_{o}$ then manifests itself as a larger $a_{j,1}$ from (5.9). This can be seen in figure 15(c) where $\bar {\sigma }_{1}$ is computed for the LNSE where only $\bar {V}$ and $\bar {V}_{y}$ are neglected. By including $\bar {U}_{x}$, $\bar {\sigma }_{1} >\tilde {\sigma }_{1}$. Since $\mathcal {L}_{SQ}^{-1}$ also amplifies $\hat {f}_{\omega }$, the magnitude of $a_{j,0}$ is also increased when $\bar {U}_{x}$ is included.

Figure 15. Comparison of $\tilde {\sigma }_{j,SQ}$ (blue) and $\sigma _{j,SQ}$ (red) (a). Here $\sigma _{j,(v\rightarrow u,w)}$ is computed for the full LNSE (blue) and $\tilde {\sigma }_{j,(v\rightarrow u,w)}$ by neglecting $\bar {U}_{x}$, $\bar {V}_{x}$, and $\bar {V}$ (red) (b). The lines of best fit are fitted using (5.8). Plots of $\sigma _{1}$ (using full LNSE), $\tilde {\sigma }_{1}$ (neglecting $\bar {V}$, $\bar {V}_{x}$, and $\bar {U}_{x}$) and $\bar {\sigma }_{1}$ (neglecting $\bar {V}$ and $\bar {V}_{x}$) in black, red and blue, respectively. The lines of best fit are fitted for $\langle \beta _{\infty } \rangle > 0$. $j = 1,$ $2$, $3$, $4$ and $5$ correspond to the circle, star, square, triangle and $\times$. These are all computed for the representative large-scale mode from figure 8.

The inclusion of the $\bar {V}$ terms provide wall-normal advection that is neither present in the ZPG TBL, nor in LNSE with artificially neglected $\bar {V}$. This can be seen in figure 16(a,b), where $| \psi _{1,u} |$ is compared using the full LNSE and the LNSE with $\bar {V}$ and $\bar {V}_{x}$ artificially set to $0$. It can be seen that the modes computed with the full LNSE reach their peak amplitude farther from the wall than the artificial modes. As the APG strength increases, the disparity between the two sets of modes increases, as can be seen by comparing the m13 and b1n modes. As expected, the ZPG modes are almost indistinguishable. Because the modes computed with stronger APG are farther from the wall, they experience a stronger streamwise varying $\bar {U}$. As a result, the modes computed with the full LNSE experience a stronger convective non-normality and thus increased amplification from the $\mathcal {L}_{SQ}^{-1}$ term. As a result, for $\sigma _{1}$, both $a_{j,1}$ and $| a_{j,0} |$ are the largest in figure 15(c), indicating a further increase in amplification when none of the non-parallel terms are neglected.

Figure 16. Plots of $| \psi _{1,u} |$ at $x_{r_{2}}$ for the representative large-scale mode. The dotted lines are computed with $\bar {V}$ and $\bar {V}_{x}$ set to $0$ while the solid lines are computed with the full LNSE. Panels (a) and (b) are separated for visibility. The lines are colour coded according to table 1.