2.1. Wall layer

The immediate response of the flow to the surface roughness is in a viscous WL with an $\mathit{O}(R^{-1/3})$ width, which is derived by balancing the inertia and viscous terms in the chordwise momentum equations in (2.3). We thus introduce in this layer a local wall-normal coordinate,

(2.18) \begin{equation} Y = R^{1/3} y = \mathit{O} (1), \end{equation}

in terms of which the wavy wall is specified as $Y_{w} = h F(x,z)$ . Note that with $h = \mathit{O} (1)$ , the chordwise and spanwise velocities of the induced disturbance are comparable to those of the background flow, both being of $\mathit{O} (R^{-1/3})$ . It follows that the total flow field in the WL is scaled as

(2.19) \begin{equation} (u,v,w,p ) = \bigl (R^{-1/3} U_{w}, R^{-2/3} V_{w} , R^{-1/3} W_{w}, R^{-2/3} P_{w} \bigr ). \end{equation}

Substitution of (2.19) into (2.3) gives boundary-layer equations,

(2.20a) \begin{align} U_{w, x} + V_{w, Y} + W_{w, z} = & 0, \end{align}

(2.20b) \begin{align} U_{w} U_{w, x} + V_{w} U_{w, Y} + W_{w} U_{w, z} = & - P_{w, x} + U_{w, Y Y}, \end{align}

(2.20c) \begin{align} 0 = & P_{w, Y}, \end{align}

(2.20d) \begin{align} U_{w} W_{w, x} + V_{w} W_{w, Y} + W_{w} W_{w, z} = & - P_{w, z} + W_{w, Y Y}. \end{align}

The no-slip and impermeability boundary conditions at the roughness surface read

(2.21a–c) \begin{align} U_{w} = V_{w} = W_{w} = 0 \quad \textrm{at} \, \, \, Y = Y_{w}; \end{align}

while the matching conditions read

(2.22a,b) \begin{align} ( U_{w}, W_{w} ) \to ( \lambda _{1}, \lambda _{3} ) Y + h ( \bar {U}_{s}, \bar {W}_{s}) \quad \textrm{as} \, \, \, Y \to \infty , \end{align}

where $\lambda _{1}$ and $\lambda _{3}$ denote the chordwise and spanwise wall shears, respectively; $\bar {U}_{s}$ and $\bar {W}_{s}$ are the chordwise and spanwise streaming caused by nonlinear interactions in the WL and are independent of $Y$ , as is shown in the supplementary material (§ S1).

In order to calculate the flow field in a rectangular domain, we apply the Prandtl transformation,

(2.23a,b) \begin{align} \bar {Y} = Y - Y_{w}, \quad \bar {V}_w = V_w - h \left ( \bar {U}_w F_x + \bar {W}_w F_z \right )\!, \end{align}

where an overbar signifies variables under the transformed coordinates. The governing equations for the transformed variables remain the same form as (2.20), while the boundary conditions are transformed into

(2.24a–c) \begin{align} \bar {U}_w = \bar {V}_w = \bar {W}_w = 0 \quad &\textrm{at} \, \, \, \bar {Y} = 0; \end{align}

(2.24d,e) \begin{align} ( \bar {U}_{w}, \bar {W}_{w} ) \to ( \lambda _{1}, \lambda _{3} ) (\bar {Y} + h F) + h ( \bar {U}_{s}, \bar {W}_{s}) \quad &\textrm{as} \, \, \, \bar {Y} \to \infty . \ \end{align}

We then decompose the total flow field into the base flow and the disturbance,

(2.25) \begin{align} (\bar {U}_w,\bar {V}_w,\bar {W}_w,\bar {P}_w) &= (\lambda _1 (\bar {Y}+h F),- h \lambda _{\textit{eff}} F_{\zeta } (\bar {Y}+h F),\lambda _3 (\bar {Y}+ h F),P_B)\nonumber\\&\quad + h ( \bar {U}^{\dagger }, \bar {V}^{\dagger }, \bar {W}^{\dagger },\bar {P}^{\dagger }), \end{align}

and express the latter as Fourier series,

(2.26) \begin{equation} (\bar {U}^{\dagger },\bar {V}^{\dagger },\bar {W}^{\dagger },\bar {P}^{\dagger }) = \sum _{n=-\infty }^{\infty } (\bar {U}_n(\bar {Y}),\bar {V}_n (\bar {Y}) ,\bar {W}_n (\bar {Y}),\bar {P}_n (\bar {Y})) \textrm{e}^{\textrm{i} n \alpha _w \zeta }, \end{equation}

where we introduce the skewed coordinate $\zeta = x + (\beta _w / \alpha _w) z$ and the effective base-flow wall shear $\lambda _{\textit{eff}} = \lambda _1 + (\beta _w / \alpha _w) \lambda _3$ . To keep the physical quantities real, we require

(2.27) \begin{equation} \bar {q}_{-n} = (\bar {q}_n)_{cc}, \end{equation}

where $\bar {q}_n$ stands for any of the Fourier components on the right-hand side of (2.26). Substituting (2.25)–(2.26) into (2.20) gives the governing equations for the Fourier components,

(2.28a) \begin{align} \textrm{i} n \alpha _w \bar {U}_n + \bar {V}_{n, \bar {Y}} + \textrm{i} n \beta _w \bar {W}_n & = 0, \end{align}

(2.28b) \begin{align} \textrm{i} n \alpha _w \lambda _{\textit{eff}} \bar {Y} \bar {U}_n + \lambda _1 \bar {V}_n + \textrm{i} n \alpha _w \bar {P}_n - \bar {U}_{n, \bar {Y} \bar {Y}} & = - h ( \lambda _{\textit{eff}} \mathcal{R}_n^x + \mathcal{N}_n^x ), \end{align}

(2.28c) \begin{align} \bar {P}_{n, \bar {Y}} & = 0, \end{align}

(2.28d) \begin{align} \textrm{i} n \alpha _w \lambda _{\textit{eff}} \bar {Y} \bar {W}_n + \lambda _3 \bar {V}_n + \textrm{i} n \beta _w \bar {P}_n - \bar {W}_{n, \bar {Y} \bar {Y}} & = - h ( \lambda _{\textit{eff}} \mathcal{R}_n^z + \mathcal{N}_n^z ), \end{align}

where

(2.29a) \begin{align} \mathcal{R}_n^x &= \left \{ \widehat {F \bar {U}_{\zeta }^{\dagger }} \right \}_n + \frac {\lambda _1}{\lambda _{\textit{eff}}} \left ( \left \{ \widehat { F_{\zeta } \bar {U}^{\dagger }} \right \}_n + \frac {\beta _w}{\alpha _w} \left \{ \widehat { F_{\zeta } \bar {W}^{\dagger }} \right \}_n \right ) - \bar {Y} \left \{ \widehat { F_{\zeta } \bar {U}_{\bar {Y}}^{\dagger } } \right \}_n - h \left \{ \widehat {F F_{\zeta } \bar {U}_{\bar {Y}}^{\dagger }} \right \}_n , \end{align}

(2.29b) \begin{align} \mathcal{R}_n^z &= \left \{ \widehat {F \bar {W}_{\zeta }^{\dagger }} \right \}_n + \frac {\lambda _3}{\lambda _{\textit{eff}}} \left ( \left \{ \widehat { F_{\zeta } \bar {U}^{\dagger }} \right \}_n + \frac {\beta _w}{\alpha _w} \left \{ \widehat { F_{\zeta } \bar {W}^{\dagger }} \right \}_n \right ) - \bar {Y} \left \{ \widehat { F_{\zeta } \bar {W}_{\bar {Y}}^{\dagger } } \right \}_n - h \left \{ \widehat {F F_{\zeta } \bar {W}_{\bar {Y}}^{\dagger }} \right \}_n , \end{align}

(2.29c) \begin{align} \mathcal{N}_n^x &= \textrm{i} n \alpha _w \left \{ \widehat {\bar {U}^{\dagger } \bar {U}^{\dagger }} \right \}_n + \left \{ \widehat {\bar {V}^{\dagger } \bar {U}^{\dagger } } \right \}_{n,\bar {Y}} + \textrm{i} n \beta _w \left \{ \widehat { \bar {W}^{\dagger } \bar {U}^{\dagger } } \right \}_n, \end{align}

(2.29d) \begin{align} \mathcal{N}_n^z &= \textrm{i} n \alpha _w \left \{ \widehat {\bar {U}^{\dagger } \bar {W}^{\dagger } } \right \}_n + \left \{ \widehat { \bar {V}^{\dagger } \bar {W}^{\dagger }} \right \}_{n,\bar {Y}} + \textrm{i} n \beta _w \left \{ \widehat {\bar {W}^{\dagger } \bar {W}^{\dagger }} \right \}_n, \end{align}

with $\widehat {\qquad }$ denoting the Fourier transform with respect to $\zeta$ and the subscript $n$ representing the $n$ th harmonic. The boundary conditions are

(2.30a–c) \begin{align} \bar {U}_{n} = - \lambda _1 F_n, \quad \bar {V}_{n} = h \lambda _{\textit{eff}} \left \{ \widehat {F F_{\zeta }} \right \}_n, \quad \bar {W}_{n} = - \lambda _3 F_n \quad \textrm{at} \, \, \, \bar {Y} = 0. \end{align}

Note that $\bar {U}_{0} (0) = \bar {W}_{0} (0) = 0$ since $F_{0} = 0$ , and $\bar {V}_{0} (0) = 0$ since $ \{ \widehat {F F_{\zeta }} \}_{n} = \textrm{i} n \alpha _{w} \{ \widehat {F^{2}} \}_{n} / 2$ . For all $n \ne 0$ Fourier components, we have the matching conditions

(2.31a,b) \begin{align} \bar {U}_{n} \to 0, \quad \bar {W}_{n} \to 0 \quad \textrm{as} \, \, \, \bar {Y} \to \infty . \end{align}

Each Fourier component of the pressure $\bar {P}_{n}$ is a constant according to (2.28c ). Combining (2.28b ) and (2.28d ) and taking the limit $\bar {Y} \to \infty$ , we obtain the expression for a $n \ne 0$ Fourier component of the pressure

(2.32) \begin{equation} \bar {P}_{n} = - \alpha _w^2 k_w^{-2} \lambda _{\textit{eff}} \bigl [ \bar {V}_{n,\infty } / (\textrm{i} n \alpha _w) + h F_n \bar {U}_{s \textit{eff}} \bigr ], \end{equation}

where $\bar {V}_{n,\infty } \equiv \bar {V}_{n} \rvert _{\bar {Y} \to \infty }$ and

(2.33a,b) \begin{align} k_{w} = \sqrt {\alpha _w^2 + \beta _w^2}, \quad \bar {U}_{s \textit{eff}} = \bar {U}_{s} + (\beta _{w} / \alpha _{w}) \bar {W}_{s}. \end{align}

The Fourier component with $n = 0$ , also referred to as the mean-flow distortion, is caused by nonlinear interactions of the roughness-induced disturbance within the WL. It exhibits a different asymptotic behaviour as $\bar {Y} \to \infty$ so that a separate formulation is required. Taking $n = 0$ in (2.28a ) gives $\bar {V}_{0, \bar {Y}} = 0$ , integration of which with respect to $\bar {Y}$ with use of the boundary condition (2.30b ) leads to $\bar {V}_{0} \equiv 0$ . Setting $n = 0$ in (2.28b ) and (2.28d ), we obtain the governing equations for $\bar {U}_{0}$ and $\bar {W}_{0}$ ,

(2.34a,b) \begin{align} \bar {U}_{0, \bar {Y} \bar {Y}} = h \bigl ( \mathcal{R}_0^x + \mathcal{N}_0^x \bigr ), \quad \bar {W}_{0, \bar {Y} \bar {Y}} = h \bigl ( \mathcal{R}_0^z + \mathcal{N}_0^z \bigr ). \end{align}

Integrating (2.34) twice, we obtain the solutions for $\bar {U}_{0}$ and $\bar {W}_{0}$ ,

(2.35a) \begin{align} \bar {U}_{0} = h \int _{0}^{\bar {Y}} ( \bar {Y} - Y_1 ) \bigl ( \mathcal{R}_0^x (Y_1) + \mathcal{N}_0^x(Y_1) \bigr ) \textrm{d} Y_1 + \bar {U}_{0,\bar {Y}} (0) \bar {Y}, \end{align}

(2.35b) \begin{align} \bar {W}_{0} = h \int _{0}^{\bar {Y}} ( \bar {Y} - Y_1 ) \bigl ( \mathcal{R}_0^z (Y_1) + \mathcal{N}_0^z(Y_1) \bigr ) \textrm{d} Y_1 + \bar {W}_{0,\bar {Y}} (0) \bar {Y}, \end{align}

where

(2.36a,b) \begin{align} \bar {U}_{0,\bar {Y}} (0) = - h \frac {\alpha _w^2}{k_w^2} \lambda _{\textit{eff}} \sum _{\substack {k=-\infty \\ k \ne 0}}^{\infty } F_{-k} \bar {V}_{k, \infty } , \ \bar {W}_{0,\bar {Y}} (0) = - h \frac {\alpha _w \beta _w}{k_w^2} \lambda _{\textit{eff}} \sum _{\substack {k=-\infty \\ k \ne 0}}^{\infty } F_{-k} \bar {V}_{k, \infty }, \end{align}

which are fixed by considering the large- $\bar {Y}$ asymptotes of $\bar {U}_{0}$ and $\bar {W}_{0}$ , as discussed in the supplementary material (§ S1). It can be shown that

(2.37a,b) \begin{align} \bar {U}_{0} \to \bar {U}_s, \quad \bar {W}_{0} \to \bar {W}_s \quad \textrm{as} \, \, \, \bar {Y} \to \infty , \end{align}

with both $\bar {U}_{s}$ and $\bar {W}_{s}$ being constants. This nonlinear streaming effect induces slip velocities which act at the bottom of the main layer, driving a mean-flow distortion of $\mathit{O} (R^{-1/3})$ there.

The WL generates a ‘blowing velocity’ of $\mathit{O} (R^{-2/3})$ as $\bar {Y} \to \infty$ , which can be derived by integrating (2.20a ) and substituting into it the matching conditions (2.24d , e ),

(2.38) \begin{equation} \bar {V}_{w} \to - h ( \lambda _{1} F_{x} + \lambda _{3} F_{z} ) \bar {Y} + \bar {V}_{t} (x, z) \quad \textrm{as} \, \, \, \bar {Y} \to \infty , \end{equation}

where $\bar {V}_{t}$ has the expression

(2.39) \begin{equation} \bar {V}_{t} = - h \sum _{n = - \infty }^{\infty } \biggl [ \int _{0}^{\infty } ( \textrm{i} n \alpha _w \bar {U}_n + \textrm{i} n \beta _w \bar {W}_n ) \textrm{d} Y \biggr ] \textrm{e}^{\textrm{i} n \alpha _w \zeta }. \end{equation}

Transforming the above quantities back to the original coordinates, we can calculate the blowing velocity of the WL,

(2.40) \begin{equation} V_w \to \bar {V}_t + h^2 \left (\lambda _{\textit{eff}} F + \bar {U}_{s \textit{eff}} \right ) F_{\zeta } \quad \textrm{as} \, \, \, Y \to \infty , \end{equation}

which acts at the bottom of the main layer to induce a response of the same order there; the first and second terms on right-hand side of (2.40) represent, respectively, the transpiration velocity induced by the viscous motion in the WL and the pure geometric deflection of the streamlines by the roughness topography. Physically, this blowing velocity can be interpreted as a blowing/suction unit on the wall, which generates three-dimensional vorticity in the main layer.

The WL dynamics is strongly nonlinear with $\mathit{O} (R^{-1/3})$ chordwise and spanwise velocities. The associated vorticity is even greater with an $\mathit{O} (1)$ amplitude, making the WL susceptible to small-scale secondary instability as will be shown in § 4.

2.2. Main layer

The main layer occupies the bulk of the boundary layer, where $y = \mathit{O} (1)$ . The $\mathit{O} (R^{-2/3})$ blowing velocity and the $\mathit{O} (R^{-1/3})$ chordwise and spanwise streaming produce corresponding responses in the main layer. The flow field there is expressed as

(2.41) \begin{align} (\boldsymbol{u}, p ) &= \big( U_{\!B}, R^{-1} V_{B}, W_{\!B},P_B \big) \nonumber \\ & \quad + R^{-1/3} \Big( U_{s}^{\langle 1 \rangle }, R^{-1} V_{s}^{\langle 1 \rangle },W_{s}^{\langle 1 \rangle }, P_{s}^{\langle 1 \rangle } \Big) + R^{-2/3} \Big( U_{s}^{\langle 2 \rangle }, R^{-1} V_{s}^{\langle 2 \rangle }, W_{s}^{\langle 2 \rangle }, P_{s}^{\langle 2 \rangle } \Big) \nonumber \\ & \quad + R^{-2/3} \sum _{\substack {n=-\infty \\ n \ne 0}}^{\infty } ( u_{n} \left ( y \right )\!, v_{n}\left ( y \right )\!, w_{n}\left ( y \right )\!, p_{n}\left ( y \right ) ) e^{i n \alpha _w \zeta } + {\cdots} , \end{align}

where the subscript $s$ denotes the main-layer streaming driven by the mean-flow distortion from the WL and the CL, which will be shown later in § 2.3. This streaming, though larger, is not our main interest since it does not affect the dynamics in the CL nor the secondary instability. We focus on the $\mathit{O} (R^{-2/3})$ forced disturbance, which is governed by

(2.42a) \begin{align} \textrm{i} n \alpha _w u_{n} + v_{n,y}+ \textrm{i} n \beta _w w_{n} = 0, \end{align}

(2.42b) \begin{align} \textrm{i} n \alpha _w U_{\!B\textit{eff}} u_{n} + v_{n} U_{\!B,y} + \textrm{i} n \alpha _w p_{n} = 0, \end{align}

(2.42c) \begin{align} \textrm{i} n \alpha _w U_{\!B\textit{eff}} v_{n} + p_{n,y} = 0, \end{align}

(2.42d) \begin{align} \textrm{i} n \alpha _w U_{\!B\textit{eff}} w_{n} + v_{n} W_{B,y} + \textrm{i} n \beta _w p_{n} = 0, \end{align}

where we have introduced the effective base flow

(2.43) \begin{equation} U_{\!B \textit{eff}} = U_B + (\beta _w / \alpha _{w}) W_B. \end{equation}

Equations in (2.42) can be reduced to the steady Rayleigh equation

(2.44) \begin{equation} v_{n, y y} - \left (n^2 k_w^2 + \frac {U_{\!B \textit{eff}, yy}}{U_{\!B \textit{eff}}} \right ) v_{n} = 0 \quad n \ne 0, \end{equation}

which is subject to the inhomogeneous boundary conditions,

(2.45a,b) \begin{align} v_{n} \to V_{w \infty ,n} \quad \textrm{as} \, \, \, y \to 0; \quad v_{n} \to 0 \quad \textrm{as} \, \, \, y \to \infty . \end{align}

Here $V_{w \infty ,n}$ represents the Fourier component on the right-hand side of (2.40).

The Rayleigh (2.44) becomes singular at the critical level $y_c$ , where $U_{\!B \textit{eff}} (y_c) = 0$ while $U_{\!B \textit{eff}, yy} (y_c) \ne 0$ in general. Near $y_c$ , a local solution to (2.44) can be constructed using the Frobenius method, which gives

(2.46) \begin{equation} v_{n} = a_{n}^{\pm } \phi _{an} (y)+b_{n} \Bigl [\phi _{bn} (y) + \frac {U_c^{\prime\prime}}{U_c^{\prime}} \ln {\left |y-y_c\right |} \phi _{an}(y) \Bigr ] \end{equation}

with

(2.47a,b) \begin{align} \phi _{an} = \tilde {y} + \frac {1}{2} \frac {U_c^{\prime \prime }}{U_c^{\prime }} \tilde {y} ^ 2 + {\cdots} , \quad \phi _{bn} = 1 + \frac {1}{2} \left (n^2 k_w^2 + \frac {U_c^{\prime \prime \prime }}{U_c^{\prime }} - 2 \frac {U_c^{\prime \prime 2}}{U_c^{\prime 2}} \right ) \tilde {y}^2 + {\cdots} , \end{align}

and $\tilde {y} = y - y_{c}$ . The subscript $c$ signifies the quantities related to $U_{\!B \textit{eff}} (y_{c})$ , and a prime denotes the derivative with respect to $y$ . There are velocity jumps $(a_n^+ - a_n^-)$ across $y_c$ , which will be determined by analysing the CL dynamics.

In the main layer, we can also express $\{ u_n,w_n,p_n \}$ in terms of $v_n$ and its derivative,

(2.48a) \begin{align} u_{n} &= \frac {\textrm{i} \alpha _w}{n k_w^2} \left (v_{n,y} - \frac {U_{\!B\textit{eff},y}}{U_{\!B\textit{eff}}} v_{n} \right ) + \frac {\textrm{i} U_{\!B,y}}{n \alpha _w U_{\!B\textit{eff}}}v_{n}, \end{align}

(2.48b) \begin{align} w_{n} &= \frac {\textrm{i} \beta _w}{n k_w^2} \left (v_{n,y} - \frac {U_{\!B\textit{eff},y}}{U_{\!B\textit{eff}}} v_{n} \right ) + \frac {\textrm{i} W_{B,y}}{n \alpha _w U_{\!B\textit{eff}}}v_{n}, \end{align}

(2.48c) \begin{align} p_{n} &= - \frac {\textrm{i} \alpha _w}{n k_w^2} \left ( v_{n,y} - \frac {U_{\!B\textit{eff},y}}{U_{\!B \textit{eff}}} v_{n} \right ) U_{\!B \textit{eff}}. \end{align}

2.3. Critical layer

To resolve the singularity of (2.44) at the critical level $y_c$ , we need to consider a CL around $y_c$ with an $\mathit{O} (R^{-1/3})$ width, within which the viscous effect is reinstated. We introduce a new local coordinate in the CL,

(2.49) \begin{equation} \tilde {Y} = R^{1/3} (y - y_c) = \mathit{O} (1). \end{equation}

In order to determine the form of expansion for the disturbance in this layer, the main-layer solutions near the critical level, (2.46) and (2.48), are rewritten in terms of $\tilde {Y}$ ,

(2.50a) \begin{align} u &\to U_{\!B} (y_c) + R^{-\frac {1}{3}} \biggl [ U_{\!B}^{\prime} (y_c) \tilde {Y} + U_{s}^{\langle 1 \rangle } (y_c) + \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {\textrm{1}}{\textrm{i} n \alpha _w} \left (\frac {\alpha _w^2}{k_w^2} - \frac {U_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \right ) b_n \frac {1}{\tilde {Y}} \textrm{e}^{\textrm{i} n \alpha _w \zeta } \biggr ] \nonumber \\ &\quad + R^{-\frac {2}{3}} \biggl [\frac {U_{\!B}^{\prime\prime} (y_\textit{c})}{2} \tilde {Y}^2 + U_{s}^{\langle 2 \rangle } (y_c) - \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {1}{\textrm{i} n \alpha _w U_c^{\prime}} \biggl ( U_{\!B}^{\prime} (y_c) a_{n}^{\pm } + U_{\!B}^{\prime\prime} (y_c) b_n \nonumber \\ &\quad + \frac {U_c^{\prime\prime}}{2} \left (\frac {\alpha _w^2}{k_w^2} - \frac {U_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \right ) b_n + \frac {U_{\!B}^{\prime}(y_c) U_c^{\prime\prime}}{U_c^{\prime}} b_n \ln \left | \tilde {Y} \right | \biggr ) \textrm{e}^{\textrm{i} n \alpha _w \zeta } \biggr ] + {\cdots} , \end{align}

(2.50b) \begin{align} v &\to R^{-1} V_{B}(y_c) + R^{-\frac {2}{3}} \sum _{\substack {n = \infty \\ n \ne 0}}^{\infty } b_n \textrm{e}^{\textrm{i} n \alpha _w \zeta } + R^{-1} \sum _{\substack {n = \infty \\ n \ne 0}}^{\infty } \biggl ( a_{n}^{\pm } + \frac {U_c^{\prime\prime}}{U_c^{\prime}} b_n \ln \left | \tilde {Y} \right | \biggr ) \tilde {Y} \textrm{e}^{\textrm{i} n \alpha _w \zeta } + {\cdots} , \end{align}

(2.50c) \begin{align} w &\to \! W_{\!B} (y_c) + R^{-\frac {1}{3}} \biggl [\! W_{\!B}^{\prime} (y_c) \tilde {Y} + W_{s}^{\langle 1 \rangle } (y_c) + \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {\textrm{1}}{\textrm{i} n \alpha _w} \!\left (\!\frac {\alpha _w \beta _w}{k_w^2} - \frac {W_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \!\right )\! b_n \frac {1}{\tilde {Y}} \textrm{e}^{\textrm{i} n \alpha _w \zeta }\! \biggr ] \nonumber \\ &\quad + R^{-\frac {2}{3}} \biggl [\frac {W_{\!B}^{\prime\prime} (y_\textit{c})}{2} \tilde {Y}^2 + W_{s}^{\langle 2 \rangle } (y_c) - \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {1}{\textrm{i} n \alpha _w U_c^{\prime}} \biggl ( W_{\!B}^{\prime} (y_c) a_{n}^{\pm } + W_{\!B}^{\prime\prime} (y_c) b_n \nonumber \\ &\quad + \frac {U_c^{\prime\prime}}{2} \left (\frac {\alpha _w \beta _w}{k_w^2} - \frac {W_{\!B}^{\prime}(y_\textit{c})}{U_c^{\prime}} \right ) b_n + \frac {W_{\!B}^{\prime}(y_c) U_c^{\prime\prime}}{U_c^{\prime}} b_n \ln \left | \tilde {Y} \right | \biggr ) \textrm{e}^{\textrm{i} n \alpha _w \zeta } \biggr ] + {\cdots} , \end{align}

(2.50d) \begin{align} p &\to P_{B} (y_c) + R^{-\frac {2}{3}} \sum _{\substack {n = -\infty \\ n \ne 0}}^{\infty } \frac {\textrm{i} \alpha _w}{n k_w^2} U_c^{\prime} b_n \textrm{e}^{\textrm{i} n \alpha _w \zeta } + {\cdots} , \end{align}

where the terms associated with $\ln R$ are excluded as they are the natural continuations of the outer solutions and remain regular in the CL. We note that $u$ and $w$ have a simple-pole singularity, while $v$ has a logarithmic singularity. Based upon the asymptotic matching principle, we express the flow field in the CL as

(2.51) \begin{align} (\boldsymbol{u}, p ) &= (U_B (y_{c}),R^{-1} V_B (y_{c}), W_B (y_{c}),P_B ) \nonumber \\ & \quad + R^{-\frac {1}{3}} ( U_{\!B}^{\prime} (y_{c}), 0, W_{\!B}^{\prime} (y_c), 0 ) \tilde {Y} + \frac {1}{2} R^{-\frac {2}{3}} ( U_{\!B}^{\prime\prime} (y_{c}), 0, W_{\!B}^{\prime\prime} (y_c), 0 ) \tilde {Y}^{2} \nonumber \\ & \quad + R^{-1/3} ( \tilde {U}, R^{-1/3} \tilde {V}, \tilde {W}, R^{-1/3} \tilde {P} ), \end{align}

where

(2.52) \begin{align} \left( \tilde {U}, R^{-1/3} \tilde {V}, \tilde {W}, R^{-1/3} \tilde {P} \right) &= \sum _{n = -\infty }^{\infty } \left [ \left(\tilde {u}_{n}^{\langle 1 \rangle } (\tilde {Y}), R^{-1/3} b_{n}, \tilde {w}_{n}^{\langle 1 \rangle } (\tilde {Y}), R^{-1/3} \tilde {p}_{n}^{\langle 1 \rangle } \right)\right . \nonumber \\ &\quad \left . + R^{-1/3} \left(\tilde {u}_{n}^{\langle 2 \rangle } (\tilde {Y}), R^{-1/3} \tilde {v}_{n}^{\langle 2 \rangle } (\tilde {Y}), \tilde {w}_{n}^{\langle 2 \rangle } (\tilde {Y}), R^{-1/3} \tilde {p}_{n}^{\langle 2 \rangle } \right) \right ]\nonumber \\ &\quad \times \textrm{e}^{\textrm{i} n \alpha _w \zeta }. \end{align}

By matching with the main-layer solution (2.50), we derive

(2.53a,b) \begin{align} b_{0} &= 0, \quad \tilde {v}_{0}^{\langle 1 \rangle } = 0; \end{align}

(2.53c–e) \begin{align} \tilde {p}_{0}^{\langle 1 \rangle } &= 0, \quad \ \tilde {p}_{n}^{\langle 1 \rangle } = (\textrm{i} \alpha _w) / (n k_w^2) U_c^{\prime} b_n \ (n \ne 0); \qquad \tilde {p}_{n}^{\langle 2 \rangle } \equiv 0. \end{align}

Substituting (2.51)–(2.52) into the N–S (2.3), we obtain, at leading order, the governing equations for each Fourier component $(n \ne 0)$ of the regularised disturbance,

(2.54a) \begin{align} \alpha _w \tilde {u}_n^{\langle 1 \rangle } + \beta _w \tilde {w}_n^{\langle 1 \rangle } &= 0, \end{align}

(2.54b) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {u}_n^{\langle 1 \rangle } + b_{n} U_{\!B}^{\prime} (y_c) + \textrm{i} n \alpha _{w} \tilde {p}_{n}^{\langle 1 \rangle } - \tilde {u}_{n,\tilde {Y} \tilde {Y}}^{\langle 1 \rangle } &= -\tilde {\mathcal{N}}_{n}^{ \left \langle 1 x \right \rangle }, \end{align}

(2.54c) \begin{align} \tilde {p}_{n,\tilde {Y}}^{\langle 1 \rangle } &= 0, \end{align}

(2.54d) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {w}_n^{\langle 1 \rangle } + b_{n} W_{B}^{\prime} (y_c) + \textrm{i} n \beta _{w} \tilde {p}_{n}^{\langle 1 \rangle } - \tilde {w}_{n,\tilde {Y} \tilde {Y}}^{\langle 1 \rangle } &= -\tilde {\mathcal{N}}_{n}^{ \left \langle 1 z \right \rangle }, \end{align}

where

(2.55a) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 1 x \right \rangle } = & \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty }b_{n-k} \tilde {u}_{k,\tilde {Y}}^{\langle 1 \rangle } + b_n \tilde {u}_{0,\tilde {Y}}^{\langle 1 \rangle } + \left (\tilde {u}_0^{\langle 1 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_0^{\langle 1 \rangle } \right ) \textrm{i} n \alpha _w \tilde {u}_n^{\langle 1 \rangle }, \end{align}

(2.55b) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 1 z \right \rangle } = & \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty }b_{n-k} \tilde {w}_{k,\tilde {Y}}^{\langle 1 \rangle } + b_n \tilde {w}_{0,\tilde {Y}}^{\langle 1 \rangle } + \left (\tilde {u}_0^{\langle 1 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_0^{\langle 1 \rangle } \right ) \textrm{i} n \alpha _w \tilde {w}_n^{\langle 1 \rangle }. \end{align}

The spanwise momentum (2.54d ) is proportional to the chordwise momentum (2.54b ) according to the continuity (2.54a ). The solution for $\tilde {u}_n^{\langle 1 \rangle } (n \ne 0 )$ is determined by solving (2.54b ) numerically subject to the boundary conditions

(2.56) \begin{equation} \tilde {u}_n^{\langle 1 \rangle } \to \textrm{i} / (n \alpha _w) \left( U_{\!B}^{\prime}(y_c) / U_{c}^{\prime} - \alpha _w^2 k_{w}^{-2} \right) b_n / \tilde {Y} \qquad \textrm{as} \, \, \, \tilde {Y} \to \pm \infty , \end{equation}

which is derived by matching the CL solution with its counterpart in the main layer.

The governing equations for the leading-order mean-flow distortion $\tilde {u}_{0}^{\langle 1 \rangle } (\tilde {Y})$ and $\tilde {w}_{0}^{\langle 1 \rangle } (\tilde {Y})$ in the CL are

(2.57a,b) \begin{align} \tilde {u}_{0, \tilde {Y} \tilde {Y}}^{\langle 1 \rangle } =\sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k, \tilde {Y}}^{\langle 1 \rangle }, \quad \tilde {w}_{0, \tilde {Y} \tilde {Y}}^{\langle 1 \rangle } = \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {w}_{k, \tilde {Y}}^{\langle 1 \rangle }. \end{align}

Integrating (2.57a ) from $-\infty$ to $\tilde {Y}$ gives

(2.58) \begin{equation} \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k}^{\langle 1 \rangle } (\tilde {Y}) + \left . \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty }, \end{equation}

and then integrating again from $0$ to $\tilde {Y}$ leads to

(2.59) \begin{equation} \tilde {u}_{0}^{\langle 1 \rangle } (\tilde {Y}) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \int _{0}^{\tilde {Y}} \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1 + \left . \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty } \tilde {Y} + \left . \tilde {u}_{0}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} = 0}. \end{equation}

The first integral term in (2.59) approaches a constant as $\tilde {Y} \to \infty$ , which can be seen by substituting into it the matching condition (2.56). To prevent the non-physical $\mathit{O} ( \tilde {Y} )$ growth as $\tilde {Y} \to \infty$ in (2.59), we need to set

(2.60) \begin{equation} \left . \tilde {u}_{0, \tilde {Y}}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty } = 0. \end{equation}

It follows from (2.59) that there is a velocity jump in the leading-order mean-flow distortion

(2.61) \begin{equation} \left . \tilde {u}_{0}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to \infty } - \left . \tilde {u}_{0}^{\langle 1 \rangle } \right \rvert _{\tilde {Y} \to -\infty } = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1, \end{equation}

where the integral must be evaluated as a Cauchy principal value. Through matching the solutions in the CL and main layer, we have the streaming jump across the CL,

(2.62) \begin{equation} U_{s}^{\langle 1 \rangle } (y_c^+) - U_{s}^{\langle 1 \rangle } (y_c^-) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1. \end{equation}

Similarly, we can derive

(2.63) \begin{equation} \tilde {w}_{0, \tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) = \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {w}_{k}^{\langle 1 \rangle } (\tilde {Y}) \end{equation}

and the spanwise streaming jump. It follows from (2.54a ), (2.58) and (2.63) that

(2.64) \begin{equation} \tilde {u}_{0,\tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) + (\beta _{w} / \alpha _{w}) \tilde {w}_{0,\tilde {Y}}^{\langle 1 \rangle } (\tilde {Y}) = 0, \end{equation}

indicating that $\tilde {u}_{0}^{\langle 1 \rangle } + (\beta _w / \alpha _w) \tilde {w}_{0}^{\langle 1 \rangle }$ should be a constant, which we denote as

(2.65) \begin{equation} U_{s \textit{cef}}^{\langle 1 \rangle } \equiv \tilde {u}_{0}^{\langle 1 \rangle } + (\beta _{w} / \alpha _{w}) \tilde {w}_{0}^{\langle 1 \rangle }. \end{equation}

With the use of (2.65) in (2.55a ), the latter is expressed as

(2.66) \begin{equation} \tilde {\mathcal{N}}_{n}^{ \left \langle 1 x \right \rangle } = \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty }b_{n-k} \tilde {u}_{k,\tilde {Y}}^{\langle 1 \rangle } + b_n \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k}^{\langle 1 \rangle } + U_{s \textit{cef}}^{\langle 1 \rangle } \textrm{i} n \alpha _w \tilde {u}_n^{\langle 1 \rangle }. \end{equation}

The governing equations for each Fourier component $(n \ne 0)$ of the second-order terms in expansion (2.52) are found to be

(2.67a) \begin{align} \textrm{i} n \alpha _w \tilde {u}_n^{\langle 2 \rangle } + \tilde {v}_{n, \tilde {Y}}^{\langle 2 \rangle } + \textrm{i} n \beta _w \tilde {w}_n^{\langle 2 \rangle } = & 0, \end{align}

(2.67b) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {u}_n^{\langle 2 \rangle } + \frac {1}{2} \textrm{i} n \alpha _w U_{c}^{\prime\prime} \tilde {Y}^2 \tilde {u}_n^{\langle 1 \rangle } + \tilde {v}_n^{\langle 2 \rangle } U_{\!B}^{\prime}(y_c) \qquad \qquad & \nonumber \\ + V_B(y_c) \tilde {u}_{n,\tilde {Y}}^{\langle 1 \rangle } + b_{n} U_{\!B}^{\prime\prime}(y_c) \tilde {Y} - \tilde {u}_{n,\tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & - \tilde {\mathcal{N}}_{n}^{ \left \langle 2 x \right \rangle }, \end{align}

(2.67c) \begin{align} \textrm{i} n \alpha _w U_c^{\prime} \tilde {Y} \tilde {w}_n^{\langle 2 \rangle } + \frac {1}{2} \textrm{i} n \alpha _w U_{c}^{\prime\prime} \tilde {Y}^2 \tilde {w}_n^{\langle 1 \rangle } + \tilde {v}_n^{\langle 2 \rangle } W_{B}^{\prime}(y_c) \qquad \qquad & \nonumber \\ + V_B(y_c) \tilde {w}_{n,\tilde {Y}}^{\langle 1 \rangle } + b_{n} W_{B}^{\prime\prime}(y_c) \tilde {Y} - \tilde {w}_{n,\tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & - \tilde {\mathcal{N}}_{n}^{ \left \langle 2 z \right \rangle }, \end{align}

where use has been made of (2.53e ) and

(2.68a) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 2 x \right \rangle } &= \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty } \left [\left (\tilde {u}_{n-k}^{\langle 2 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_{n-k}^{\langle 2 \rangle } \right ) \textrm{i} k \alpha _w \tilde {u}_{k}^{\langle 1 \rangle } + b_{n-k} \tilde {u}_{k,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_{n-k}^{\langle 2 \rangle } \tilde {u}_{k,\tilde {Y}}^{\langle 1 \rangle } \right ] + U_{s \textit{cef}}^{\langle 1 \rangle } \textrm{i} n \alpha _w \tilde {u}_{n}^{ \langle 2 \rangle } \nonumber \\ &\quad + \tilde {u}_{0\textit{eff}}^{\langle 2 \rangle } \textrm{i} n \alpha _w \tilde {u}_{n}^{ \langle 1 \rangle } + b_n \tilde {u}_{0,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_n^{\langle 2 \rangle } \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{-k} \tilde {u}_{k}^{\langle 1 \rangle }, \end{align}

(2.68b) \begin{align} \tilde {\mathcal{N}}_{n}^{ \left \langle 2 z \right \rangle } &= \sum _{\substack {k=-\infty \\ k \ne 0, n }}^{\infty } \left [\left (\tilde {u}_{n-k}^{\langle 2 \rangle } + \frac {\beta _w}{\alpha _w} \tilde {w}_{n-k}^{\langle 2 \rangle } \right ) \textrm{i} k \alpha _w \tilde {w}_{k}^{\langle 1 \rangle } + b_{n-k} \tilde {w}_{k,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_{n-k}^{\langle 2 \rangle } \tilde {w}_{k,\tilde {Y}}^{\langle 1 \rangle } \right ] + U_{s \textit{cef}}^{\langle 1 \rangle } \textrm{i} n \alpha _w \tilde {w}_{n}^{ \langle 2 \rangle } \nonumber \\ &\quad + \tilde {u}_{0\textit{eff}}^{\langle 2 \rangle } \textrm{i} n \alpha _w \tilde {w}_{n}^{ \langle 1 \rangle } + b_n \tilde {w}_{0,\tilde {Y}}^{\langle 2 \rangle } + \tilde {v}_n^{\langle 2 \rangle } \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } b_{-k} \tilde {w}_{k}^{\langle 1 \rangle } \end{align}

with $\tilde {u}_{0\textit{eff}}^{\langle 2 \rangle } = \tilde {u}_{0}^{\langle 2 \rangle } + (\beta _{w} / \alpha _{w}) \tilde {w}_{0}^{\langle 2 \rangle }$ . Combining (2.67b ) and (2.67c ), and making use of (2.67a ) followed by differentiating the resulting equation with respect to $\tilde {Y}$ , we obtain

(2.69) \begin{equation} Q_{n,\tilde {Y} \tilde {Y}} - \textrm{i} n \alpha _w U_{c}^{\prime} \left ( \tilde {Y} + \frac {U_{s \textit{cef}}^{\langle 1 \rangle }}{U_{c}^{\prime}} \right ) Q_{n} + \textrm{i} n \alpha _w U_c^{\prime\prime} b_{n} = n \sum _{\substack {k = -\infty \\ k \ne 0, n}}^{\infty } \frac {1}{k} b_{n-k} Q_{k, \tilde {Y}} - \textrm{i} n \alpha _w b_n \tilde {u}_{0\textit{eff},\tilde {Y} \tilde {Y}}, \end{equation}

where $Q_{n} = \tilde {v}_{n, \tilde {Y} \tilde {Y}}^{\langle 2 \rangle }$ .

The mean-flow distortions at $\mathit{O} (R^{-2/3})$ are governed by

(2.70a) \begin{align} \tilde {u}_{0, \tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {u}_{k}^{\langle 1 \rangle } + \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \left ( \tilde {v}_{-k}^{\langle 2 \rangle } \tilde {u}_{k}^{\langle 1 \rangle } + b_{- k} \tilde {u}_{k}^{\langle 2 \rangle } \right )_{\tilde {Y}}, \end{align}

(2.70b) \begin{align} \tilde {w}_{0, \tilde {Y} \tilde {Y}}^{\langle 2 \rangle } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} \tilde {w}_{k}^{\langle 1 \rangle } + \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \left ( \tilde {v}_{-k}^{\langle 2 \rangle } \tilde {w}_{k}^{\langle 1 \rangle } + b_{- k} \tilde {w}_{k}^{\langle 2 \rangle } \right )_{\tilde {Y}}, \end{align}

which can be integrated from $- \infty$ to $\infty$ to show that

(2.71a) \begin{align} \left . \tilde {u}_{0, \tilde {Y}}^{\langle 2 \rangle } \right \rvert _{-\infty }^{\infty } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {u}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1 + \frac {\textrm{i} \alpha _w}{k_w^2} \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \frac {1}{k} \left ( a_{k}^{+} - a_{k}^{-} \right ) b_{-k}, \end{align}

(2.71b) \begin{align} \left . \tilde {w}_{0, \tilde {Y}}^{\langle 2 \rangle } \right \rvert _{-\infty }^{\infty } = & V_{B}(y_c) \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } b_{- k} -\!\!\!\!\!\int _{-\infty }^{\infty } \tilde {w}_{k}^{\langle 1 \rangle } (Y_1) \textrm{d} Y_1 + \frac {\textrm{i} \beta _w}{k_w^2} \sum _{\substack {k = - \infty \\ k \ne 0}}^{\infty } \frac {1}{k} \left ( a_{k}^{+} - a_{k}^{-} \right ) b_{-k}, \quad \end{align}

where use has been made of (2.50). Equation (2.71) are the gradient jumps of the mean-flow distortion at $\mathit{O} (R^{-2/3})$ . We combine (2.70a ) and (2.70b ) using the continuity (2.54a ) and (2.67a ) to obtain

(2.72) \begin{equation} \tilde {u}_{0\textit{eff},\tilde {Y} \tilde {Y}} = \sum _{\substack {k = - \infty k \ne 0}}^{\infty } \frac {\textrm{i}}{k \alpha _w} b_{-k} Q_{k}, \end{equation}

which can be inserted into the right-hand side of (2.69), giving

(2.73) \begin{equation} Q_{n,\tilde {Y} \tilde {Y}} - \textrm{i} n \alpha _w U_{c}^{\prime} \left ( \tilde {Y} + \frac {U_{s \textit{eff},c }^{\langle 1 \rangle }}{U_{c}^{\prime}} \right ) Q_{n} + \textrm{i} n \alpha _w U_c^{\prime\prime} b_{n} = n \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } \frac {1}{k} \bigl ( b_{n-k} Q_{k, \tilde {Y}} + b_n b_{-k} Q_{k} \bigr ). \qquad \end{equation}

To eliminate the effect of the steady streaming $U_{s \textit{cef}}^{\langle 1 \rangle }$ from the calculations of the CL, we perform a coordinate shift in the transverse direction

(2.74) \begin{equation} \mathcal{Y} = \tilde {Y} + U_{s \textit{cef}}^{\langle 1 \rangle } / U_{c}^{\prime}, \end{equation}

and rewrite (2.73) as

(2.75) \begin{equation} Q_{n,\mathcal{Y} \mathcal{Y}} - \textrm{i} n \alpha _w U_{c}^{\prime} \mathcal{Y} Q_{n} + \textrm{i} n \alpha _w U_c^{\prime \prime } b_{n} = n \sum _{\substack {k = -\infty \\ k \ne 0}}^{\infty } \frac {1}{k} \bigl ( b_{n-k} Q_{k, \mathcal{Y}} + b_n b_{-k} Q_{k} \bigr ). \end{equation}

The above system of coupled nonlinear equations is to be solved subject to the far-field boundary conditions

(2.76) \begin{equation} Q_{n} \to \big(U_{c}^{\prime \prime } / U_{c}^{\prime }\big) b_{n} / \mathcal{Y} \quad \textrm{as} \, \, \, \mathcal{Y} \to \pm \infty . \end{equation}

We denote the velocity jumps across the CL as

(2.77) \begin{equation} J_{n} = a_{n}^{+} - a_{n}^{-} \quad (n \ne 0), \end{equation}

and matching $Q_n$ with the solution in the main layer gives

(2.78) \begin{equation} J_{n} = \,\kern-1.5pt\hbox{--}\kern0.5pt\!\!\!\!\!\int _{- \infty }^{\infty } Q_{n} \left ( Y_1 \right ) \textrm{d} Y_1. \end{equation}

The governing equations and boundary conditions for $Q_{n}$ (2.75)–(2.76) and the matching condition (2.78) along with the Rayleigh (2.44) and the Frobenius solution near the critical level (2.46) are used simultaneously to resolve the coupled dynamics between the main layer and CL.

The CL is strongly nonlinear as (2.54b ) and (2.73) indicate. All harmonics of the roughness-induced disturbance are generated at the same order. Importantly, due to three-dimensionality of the disturbance and the algebraic singularity of their chordwise and spanwise velocities, the regularised disturbance velocity in the CL has a larger amplitude of $\mathit{O} (R^{-1/3})$ than its $\mathit{O}(R^{-2/3})$ counterpart in the main layer. Furthermore, the associated vorticity in the CL is even larger with an $\mathit{O} (1)$ amplitude. This can give rise to high-frequency short-wavelength secondary instability, as will be shown in § 4.

To summarise, the WL converts surface displacement into a blowing velocity of $\mathit{O} (R^{-2/3})$ , which drives a forced disturbance in the main layer. The forced disturbance is governed by the steady Rayleigh equation with a singularity at the critical level $y = y_{c}$ , which is regularised by introducing viscous and nonlinear effects in the CL. The disturbance inside the CL has a larger, $\mathit{O} (R^{-1/3})$ , amplitude, and leads to $\mathit{O} (1)$ vorticities. Physically, the effect of such surface roughness is significant near the wall and around the critical level $y_{c}$ , where the flows are nonlinearly distorted.

3.1. Numerical methods

For the roughness-induced disturbance in the WL, we solve (2.28) using the second-order finite-difference scheme to calculate $n \ne 0$ Fourier components. Only half of the Fourier components $(n \gt 0)$ need to be computed with the other half $(n \lt 0)$ being given by (2.27). The $n = 0$ Fourier component is calculated by integrating (2.35b ) using the trapezoidal formula. The roughness interaction and nonlinear terms $\{ \mathcal{R}_{n}^{x}, \mathcal{R}_{n}^{z},\mathcal{N}_{n}^{x},\mathcal{N}_{n}^{z} \}$ on the right-hand side of (2.28) are evaluated using a pseudospectral method, which transforms the quantities in the spectral space to the physical space to carry out the multiplications and then transforms the products back to the spectral space. The $3/2$ rule (Canuto Reference Canuto1988) is applied to eliminate the aliasing error. An iteration process is needed for these nonlinear calculations. We initially simplify (2.28) by setting the roughness interaction and nonlinear terms on the right-hand side to zero, which enables us to compute numerically the linear solution or to derive the analytical solution. The linear solution then serves as the initial guess for starting the iteration process solving the nonlinear problem for a small roughness height $h$ . The iteration continues until the criterion

(3.1) \begin{equation} \| \boldsymbol{\varPhi } - \boldsymbol{\varPhi }_{\textit{old}} \|_{2} / \| \boldsymbol{\varPhi } \|_{2} \le 10^{-7} \end{equation}

is satisfied, where $\boldsymbol{\varPhi }$ is the solution matrix consisting of all Fourier components of the velocity $\{ \bar {U}^{\dagger },\bar {V}^{\dagger },\bar {W}^{\dagger } \}$ . Subsequently, we employ a parameter continuation method, incrementally increasing the roughness height and using the convergent solution as the initial guess to advance the calculation for larger $h$ . To enhance the convergence, we use the under-relaxation technique (Xu, Zhang & Wu Reference Xu, Zhang and Wu2017).

For the main-layer dynamics, our primary interest is in the forced disturbance driven by the blowing velocity (2.40) from the WL motion. The forced disturbance is governed by the steady Rayleigh equation (2.44) subject to inhomogeneous boundary conditions (2.45). Apart from the singularity at $y = y_c$ , the Rayleigh equation is also singular at the lower boundary $y = 0$ because the effective base flow satisfies the no-slip boundary condition while $U_{\!B \textit{eff}, yy} (0) \ne 0$ . We need to construct a Frobenius solution to (2.44) at $y = 0$ similar to that near $y = y_c$ , namely

(3.2) \begin{equation} v_{wn} = a_{n}^{\left \langle w \right \rangle } \phi _{an}^{\left \langle w \right \rangle } (y) + b_{n}^{\left \langle w \right \rangle } \Bigl [ \phi _{b n}^{\left \langle w \right \rangle }(y) + \frac {U_{w}^{\prime\prime}}{U_{w}^{\prime}} \ln { \left |y \right | } \phi _{an}^{\left \langle w \right \rangle }(y) \Bigr ], \end{equation}

where

(3.3a,b) \begin{gather} \phi _{an}^{\left \langle w \right \rangle } = y + \frac {1}{2} \frac {U_w^{\prime \prime }}{U_w^{\prime }} y ^ 2 + {\cdots} , \quad \phi _{bn}^{\left \langle w \right \rangle } = 1 + \frac {1}{2} \left (n^2 k_w^2 + \frac {U_w^{\prime \prime \prime }}{U_w^{\prime }} - 2 \frac {U_w^{\prime \prime 2}}{U_w^{\prime \prime 2}} \right )y^2 + {\cdots} , \end{gather}

with the subscript $w$ of $U_{w}$ signifying the quantities related to $U_{\!B \textit{eff}}$ and its derivatives evaluated at $y = 0$ . Imposing the lower boundary condition (2.45a ) on (3.2) gives

(3.4) \begin{equation} b_{n}^{\left \langle w \right \rangle } = V_{w \infty , n}. \end{equation}

We start by calculating the linear solution $v_{n}$ in the main layer. The solution to the steady Rayleigh equation (2.44) has the far-field asymptote, $v_{n} \to C_{n}^{\left \langle \infty \right \rangle } \textrm{exp}(- n k_w) y$ as $y \to \infty$ , where $C_{n}^{\left \langle \infty \right \rangle }$ is an unknown constant. We first use the MATLAB solver ODE15s to integrate (2.44) downwards from a large $y$ position, denoted as $y_{\infty }$ , to $y_{c} + \delta y$ with the initial conditions

(3.5a,b) \begin{gather} \check {v}_{n} = \exp \left (- n k_{w} y_{\infty } \right )\!, \quad \check {v}_{n}' = - n k_{w} \exp \left (- n k_{w} y_{\infty } \right )\!. \end{gather}

The resulting numerical solution at $y_{c} + \delta y$ , denoted as $\{ \check {v}_{n} (y_{c} + \delta y), \check {v}_{n}' (y_{c} + \delta y) \}$ , is equated to, after being multiplied by $C_{n}^{\left \langle \infty \right \rangle }$ , the local Frobenius solution (2.46) at $y = y_{c} + \delta y$ ,

(3.6a) \begin{align} C_{n}^{\left \langle \infty \right \rangle } \check {v}_{n} (y_{c} + \delta y) = a_{n}^{+} \phi _{an} (y_{c} + \delta y) + b_{n} \varphi _{bn} (y_{c} + \delta y), \end{align}

(3.6b) \begin{align} C_{n}^{\left \langle \infty \right \rangle } \check {v}_{n}' (y_{c} + \delta y) = a_{n}^{+} \phi _{an}^{\prime} (y_{c} + \delta y) + b_{n} \varphi _{bn}^{\prime} (y_{c} + \delta y), \end{align}

where $\varphi _{bn} (y) = \phi _{bn} (y) + (U_{c}^{\prime\prime} / U_{c}^{\prime}) \ln |y - y_{c}| \phi _{an} (y)$ . We also integrate (2.44) downwards from $y_{c} - \delta y$ to $\delta y$ with the initial conditions,

(3.7a,b) \begin{gather} \hat {v}_{an} = \phi _{an} (y_c - \delta y), \quad \hat {v}_{an}^{\prime} = \phi _{an}^{\prime} (y_c - \delta y), \end{gather}

and

(3.8a,b) \begin{gather} \hat {v}_{bn} = \varphi _{bn} (y_c - \delta y), \quad \hat {v}_{bn}^{\prime} = \varphi _{bn}^{\prime} (y_c - \delta y), \end{gather}

respectively, and the numerical results at $y = \delta y$ are denoted as $\{ \hat {v}_{an} (\delta y), \hat {v}_{an}^{\prime} (\delta y) \}$ and $\{ \hat {v}_{bn} (\delta y), \hat {v}_{bn}^{\prime} (\delta y) \}$ . The linear combination of the above two, $a_{n}^{-} \hat {v}_{an} + b_{n} \hat {v}_{bn}$ , gives the total solution. Equating it and its derivative to the local Frobenius solution (3.2) at $y = \delta y$ , we obtain

(3.9a) \begin{align} a_{n}^{-} \hat {v}_{an} (\delta y) + b_{n} \hat {v}_{bn} (\delta y) = a_{n}^{\left \langle w \right \rangle } \phi _{an}^{\left \langle w \right \rangle } (\delta y) + b_{n}^{\left \langle w \right \rangle } \varphi _{wn} (\delta y), \end{align}

(3.9b) \begin{align} a_{n}^{-} \hat {v}_{an}^{\prime} (\delta y) + b_{n} \hat {v}_{bn}^{\prime} (\delta y) = a_{n}^{\left \langle w \right \rangle } \phi _{an}^{\left \langle w \right \rangle }{'} (\delta y) + b_{n}^{\left \langle w \right \rangle } \varphi _{wn}^{\prime} (\delta y), \end{align}

where $\varphi _{wn} (y) = \phi _{bn}^{\left \langle w \right \rangle } ( y) + (U_{w}^{\prime\prime} / U_{w}^{\prime}) \ln | y | \phi _{an}^{\left \langle w \right \rangle } ( y)$ . Using (3.6) and (3.9) as well as the linear velocity jump condition across the CL,

(3.10) \begin{equation} a_{n}^{+} - a_{n}^{-} = \mathcal{K}_{c} \pi b_{n}, \end{equation}

where $\mathcal{K}_{c} = \textrm{i} (U_{c}^{\prime\prime} / U_{c}^{\prime}) \textrm{sign} (\alpha _{w} U_{c}^{\prime})$ , we can form a linear system to solve the coefficients $\{ a_n^{\pm },b_n,a_n^{\left \langle w \right \rangle } \}$ in the Frobenius solutions (2.46) and (3.2) as well as the coefficient in the large- $y$ asymptote, $C_{n}^{\left \langle \infty \right \rangle }$ , for different modes $n$ ( $n \ne 0$ ),

(3.11) \begin{equation} \boldsymbol{L}_{n} \boldsymbol{x}_{n} = \boldsymbol{B}_{n}, \end{equation}

where

(3.12a) \begin{align} \boldsymbol{L}_{n} & = \begin{pmatrix} \phi _{an}(y_c + \delta y) & \varphi _{bn} (y_{c} + \delta y) & 0 & 0 & - \check {v}_n(y_c + \delta y) \\ \phi _{an}^{\prime}(y_c + \delta y) & \varphi _{bn}^{\prime} (y_{c} + \delta y) & 0 & 0 & - \check {v}_{n}' (y_c + \delta y) \\ 1 & - \mathcal{K}_{c} \pi & -1 & 0 & 0 \\ 0 & \hat {v}_{b n} (\delta y) & \hat {v}_{an} (\delta y) & -\phi _{an}^{\left \langle w \right \rangle } (\delta y) & 0 \\ 0 & \hat {v}_{b n}' (\delta y) & \hat {v}_{an}^{\prime} (\delta y) & - \phi _{an}^{\left \langle w \right \rangle }{'} (\delta y) & 0 \end{pmatrix}\!, \qquad \end{align}

(3.12b) \begin{align} \boldsymbol{x}_{n} &= {\begin{pmatrix} a_n^+ & b_n & a_n^- & a_{n}^{\left \langle w \right \rangle } & C_{n}^{\left \langle \infty \right \rangle } \end{pmatrix}}^{\textrm{T}}, \end{align}

(3.12c) \begin{align} \boldsymbol{B}_{n} &= \begin{pmatrix} 0 & 0 & 0 & b_{n}^{\left \langle w \right \rangle } \varphi _{wn} (\delta y) & b_{n}^{\left \langle w \right \rangle } \varphi _{wn}^{\prime} (\delta y) \end{pmatrix}^{\textrm{T}}. \end{align}

After obtaining the linear main-layer solution, we can calculate the linear solution for $Q_n$ in the CL by solving (2.75) with the right-hand side terms being set to $0$ . The derivative in (2.75) is approximated using the second-order finite-difference scheme. The linear solution for $Q_{n}$ can also be found analytically, but the detail is neglected here.

The linear results for $\{ b_{n}, Q_{n} \}$ are used as the initial guesses for nonlinear calculations for a small roughness height. Similar to solving the nonlinear equations in the WL, we adopt a parameter continuation method to incrementally increase the roughness height, using the convergent result for the previous slightly smaller height as the initial guesses to start iterations. The nonlinear terms are evaluated using the pseudospectral method. After obtaining $Q_{n}$ , we can calculate the velocity jump $J_{n}$ according to (2.78). With the nonlinear jump $J_n$ , we use a modified linear system compared with (3.11) to update $\boldsymbol{x}_{n}$ ,

(3.13) \begin{equation} \boldsymbol{L}_{n} \boldsymbol{x}_{n} = \boldsymbol{B}_{n}^{\left \langle N \right \rangle }, \end{equation}

where

(3.14) \begin{equation} \boldsymbol{B}_{n}^{\left \langle N \right \rangle } = \begin{pmatrix} 0 & 0 & J_{n} - \mathcal{K}_{c} \pi b_{n,0} & b_{n}^{\left \langle w \right \rangle } \varphi _{wn} (\delta y) & b_{n}^{\left \langle w \right \rangle } \varphi _{wn}^{\prime} (\delta y) \end{pmatrix}^{\textrm{T}}, \end{equation}

and $b_{n,0}$ is the value of $b_{n}$ calculated in the previous iteration. The blowing velocity $V_{w\infty ,n}$ for the current roughness height is applied on the last two entries of (3.14) according to (3.4). We implement the iteration for the above calculations to resolve the coupled CL and main-layer dynamics until the following criteria are satisfied:

(3.15a,b) \begin{gather} \| \boldsymbol{b} - \boldsymbol{b}_{\textit{old}} \|_2 / \| \boldsymbol{b} \|_2 \le 10^{-7} \quad \textrm{and} \quad \| \boldsymbol{Q} - \boldsymbol{Q}_{\textit{old}} \|_2 / \| \boldsymbol{Q} \|_2 \le 10^{-7}, \end{gather}

where $\boldsymbol{b} = \{ b_{n} \}$ and $\boldsymbol{Q} = \{ Q_{n} \}$ are the solution vector and matrix consisting of all Fourier components, respectively. We then calculate $\tilde {u}_{n}^{\langle 1 \rangle }$ in the CL by solving numerically (2.54b ), and further obtain $\tilde {w}_{n}^{\langle 1 \rangle }$ using the continuity (2.54a ).

3.2. Numerical results for disturbances in the WL

We consider the simplest roughness configuration with $F_{1} = F_{-1} = 1$ and $F_{n} = 0 \ (n \ne \pm 1)$ in (2.2) for our numerical calculations. A diagram of the surface roughness is shown in figure 2. Two pairs of roughness wavenumbers, Case I with $(\alpha _{w},\beta _{w}) = (-0.525,0.606)$ and Case II with $(\alpha _{w},\beta _{w}) = (-0.262,0.303)$ , are chosen for numerical calculations. They were chosen based on simple guiding principles: they are more-or-less of $\mathit{O} (1)$ as assumed in the theory, and fall within the wavenumber band where stationary cross-flow vortices arise. The first pair was somewhat close to that of the neutral mode of the stationary cross-flow vortices ( $(\alpha , \beta ) = (-0.6376,0.7621)$ ), while the second took the half-values of the first. We present numerical results for Case I here, while those for Case II are shown in the supplementary material (§ S2). For the base FSC flow, we set the spanwise slip velocity $W_{\infty } = 1$ and the acceleration parameter $m = 0.5$ , as it characterises the portion of typical swept-wing boundary-layer flows where a strong favourable pressure gradient is present.

Figure 2. Diagram of the surface roughness used in calculations.

We start our numerical calculations in the WL from a small roughness height ( $h = 0.1$ ) and gradually increase it. In the wall-normal direction, the domain is $\bar {Y} \in [0,30]$ with the normal grid size being set as $\Delta \bar {Y} = 0.05$ , and 17 Fourier modes $(-8 \le n \le 8)$ are included. These discretisation choices were chosen based on resolution checks, with due consideration of managing computational cost. This adequacy was supported by several observations: the results near the truncated computational domain matched the corresponding asymptotes; the current results agree with those calculated using a halved mesh size; the magnitudes of the highest Fourier mode ( $n = \pm 8$ ) remain below the iteration stopping criteria ( $10^{-7}$ ). Our numerical codes can converge up to $h = 0.825$ for Case I roughness. Figures 3 and 4 show the profiles of several harmonics of the chordwise and wall-normal velocities of the disturbance in the WL calculated for different roughness heights. The profiles of the spanwise velocities are similar to those of the chordwise velocities, and hence are omitted here. For the fundamental mode $n = 1$ , the profiles for different roughness heights, normalised by $h$ (see (2.25)), almost overlap. This is due to the fact that the flow is primarily disturbed by the fundamental mode of the roughness shape function $F$ . For higher harmonics of the velocity components, their profiles also nearly overlap when appropriately normalised, i.e. by dividing $\bar {U}_{2}$ here by $h$ and $\bar {U}_{3}$ by $h^{2}$ . Figure 4 shows that the Fourier components of the wall-normal velocity approach different constants as $\bar {Y} \to \infty$ . These constants (or equivalently, the combination of the chordwise and spanwise velocity components through the continuity (2.28a )) are used to calculate the transpiration velocity $\bar {V}_t$ , and furthermore the blowing velocity $V_w |_{\bar {Y} \to \infty }$ , which drives the $\mathit{O} (R^{- 2/3})$ disturbances in the main layer.

Figure 3. Profiles of the chordwise velocity harmonics $\bar {U}_{n}$ in the WL for different roughness height.

Figure 4. Profiles of the wall-normal velocity harmonics $\bar {V}_{n}$ in the WL for different roughness height.

Figure 5 shows the mean-flow distortions (2.35b ) in the WL. They are caused by nonlinear interactions of the disturbance, and approach constants as $\bar {Y} \to \infty$ . Such streaming in the nonlinear WL drives the $\mathit{O} (R^{-1/3})$ mean-flow distortion in the main layer. Note that the chordwise velocity is approximately four times as large as the spanwise velocity.

Figure 5. Profiles of the mean-flow distortions in the WL for different roughness height.

The $\mathit{O} (R^{-1/3})$ disturbance velocity components, along with the base flow, correspond to $\mathit{O} (1)$ vorticities in the WL, namely

(3.16) \begin{equation} \left(\partial _{x}, R^{1/3} \partial _{\bar {Y}}, \partial _{z}\right) \times \left( R^{-1/3} \bar {U}_{w}, R^{-2/3} \bar {V}_{w}, R^{-1/3} \bar {W}_{w}\right) = \left(\bar {W}_{w, \bar {Y}},0, -\bar {U}_{w, \bar {Y}}\right) + {\cdots} . \end{equation}

To aid the study of the secondary instability of the distorted WL to be conducted in § 4, we rotate the $x$ – $z$ plane of the original coordinate system so that the $x$ -axis aligns with the direction of $(\bar {\alpha },\bar {\beta })$ , where $\{ \bar {\alpha }, \bar {\beta } \}$ , introduced in (4.16), are wavenumbers of the secondary instability mode. The chordwise and spanwise velocity components, $\{ \bar {U}_{w}, \bar {W}_{w} \}$ , are projected onto the new coordinate axes, as $\{ \bar {\mathfrak{U}}_{w}, \bar {\mathfrak{W}}_{w} \}$ ,

(3.17) \begin{equation} \begin{pmatrix} \bar {\mathfrak{U}}_{w} \\ \bar {\mathfrak{W}}_{w} \end{pmatrix} = \begin{pmatrix} \cos \bar {\theta } & \sin \bar {\theta } \\ - \sin \bar {\theta } & \cos \bar {\theta } \end{pmatrix} \begin{pmatrix} \bar {U}_{w} \\ \bar {W}_{w} \end{pmatrix} \end{equation}

with $\bar {\theta } = \arctan (\bar {\beta } / \bar {\alpha })$ . The vorticity along the direction of $(\bar {\alpha },\bar {\beta })$ can be calculated as

(3.18) \begin{equation} \bar {\varOmega } = \bar {\mathfrak{U}}_{w,\bar {Y}} + {\cdots} = - \left ( \bar {\alpha } \Big / \sqrt {\bar {\alpha }^{2} + \bar {\beta }^{2}} \right ) \Bigl [ ( \bar {\beta } / \bar {\alpha } ) \bar {U}_{w} - \bar {W}_{w} \Bigr ]_{\bar {Y}} + {\cdots} . \end{equation}

Figure 6 displays contours of this skewed vorticity $\bar {\varOmega }$ without the contributions from base-flow wall shears for two roughness heights. The overall features are similar, but the vorticity induced by higher roughness is stronger, and acquires a more or less $\mathit{O} (1)$ intensity, as expected. The vorticity in the direction perpendicular to $(\bar {\alpha },\bar {\beta })$ (the direction of $(-\bar {\beta }, \bar {\alpha })$ ) controls the secondary instability, and its contours for the roughness with $h = 0.825$ will be shown later in § 5.2.1, where instability results are presented.

Figure 6. Contours of the skewed WL vorticity $\bar {\varOmega }$ without the contributions from base-flow wall shears in Case I: $(a)$ $h = 0.2$ and $(b)$ $h = 0.825$ .

The blowing velocity from the WL is of importance as it drives the response in the main layer and CL. Figure 7 shows the modulus of the first three Fourier components in the blowing velocity versus the roughness height. The fundamental mode has a much more significant modulus than higher harmonics.

Figure 7. Modulus of Fourier components in the blowing velocity versus the roughness height.

3.3. Numerical results for disturbances in the main and critical layers

For the calculations pertaining to the main layer, we first need to ascertain the position of the critical level $y_{c}$ , which is $y_{c} = 2.22081$ . Then we calculate the linear solution to the Rayleigh equation (2.44). In practice, we solve (2.44) subject to unity blowing velocity, $v_{n} (0) = b_{n}^{\left \langle w \right \rangle } = 1$ . Figure 8 presents the profiles of the first three harmonics of $v_{n} / V_{w \infty , n}$ , the linear solutions normalised by the blowing velocity. The resulting solution, multiplied by corresponding Fourier components of the blowing velocity induced by the roughness with a small height, serves as initial guesses to start the nonlinear calculations for the coupled main-layer and CL dynamics. Numerically, the linear solutions in the main layer are obtained by solving the system (3.11).

Figure 8. The normalised solutions to Rayleigh equation (2.44).

Next, we consider the coupled main-layer and CL dynamics by solving the system (3.13) along with (2.75)–(2.78). The coupling is facilitated through the velocity jumps $J_{n}$ across the CL, as shown in (3.14). Since the mean-flow distortion in the main layer is not our primary concern, and its effect on the CL dynamics amounts to a simple coordinate shift (2.74), we set $U_{s \textit{cef}}^{\langle 1 \rangle } = 0$ in our calculations and revert the wall-normal coordinate from $\mathcal{Y}$ to $\tilde {Y}$ in the CL. The system governing CL dynamics is truncated in the region $-30 \le \tilde {Y} \le 30$ , with the mesh size $\Delta \tilde {Y} = 0.05$ . The Fourier series for $b_{n}$ and $Q_{n}$ consist of 17 harmonics, i.e. $-8 \le N \le 8$ . These numerical settings were determined with both accuracy and computational efficiency in mind. Our numerical codes for this coupled dynamics can converge up to $h = 0.43$ for Case I roughness. As with the WL calculations, the nonlinear numerical codes fail to converge beyond a certain threshold, possibly due to the limitation of the pseudospectrum-based iteration; for higher $h$ , full Newton iteration may be needed. The dynamics within the CL changes from nearly linear to strongly nonlinear as the height of surface roughness is increased.

In our calculations, the Fourier components of the roughness shape function are only non-zero for $n = \pm 1$ , and the wavenumbers of Case I roughness are fairly close to those of the neutral mode of the stationary cross-flow vortices, which results in particularly significant response in the $n = \pm 1$ Fourier components of the disturbance. The nonlinear interactions of the fundamental modes generate the second and higher harmonics. We monitor only the first two harmonics of the velocity jumps $J_{n}$ . To scale out the coefficient $\mathcal{K}_{c}$ and $b_{n}$ in (3.10), we introduce a normalised jump,

(3.19) \begin{equation} \tilde {J}_{n} = J_{n} / (\mathcal{K}_{c} b_{n}), \quad n = 1, 2. \end{equation}

Figure 9 illustrates the variations of $\tilde {J}_1$ and $\tilde {J}_{2}$ versus the roughness height $h$ . When $h$ is small, $\tilde {J}_1 \to \pi$ ; whereas for sufficiently large $h$ , $\textrm{Re} \{ \tilde {J}_1 \}$ approaches zero but $\textrm{Im} \{ \tilde {J}_1 \}$ does not, which indicates a changeover from a linear viscous CL to a nonlinear viscous CL. As shown in the zoomed-in portion of figure 9(a), $\tilde {J}_{2}$ approaches a constant near $\pi$ , but not exactly $\pi$ , when $h \to 0$ due to nonlinearity. The variation of the velocity jump with $h$ can be non-monotonic.

Figure 9. Normalised velocity jumps $\tilde {J}_{1}$ and $\tilde {J}_{2}$ across the CL.

Figure 10 shows the profiles of different harmonics of $Q_{n}$ in the CL. The fundamental and second harmonic have comparable amplitudes, while the third harmonic is appreciably weaker. Higher harmonics appear to be more oscillatory since it is more substantially affected by nonlinearity. After obtaining $Q_{n}$ and $b_{n}$ , we can evaluate the nonlinear terms on the right-hand side of (2.54b ) and calculate $\tilde {u}_n^{\langle 1 \rangle }$ $(n \ne 0)$ , the profiles of which are shown in figure 11. The spanwise velocity $\tilde {w}_n^{\langle 1 \rangle }$ $(n \ne 0)$ is proportional to $\tilde {u}_{n}^{\langle 1 \rangle }$ as indicated by (2.54a ). The fundamental has a much greater amplitude than the second and third harmonics. This is due to the fact that its wavenumbers are close to those of the neutral stationary cross-flow mode, which leads to a near resonance. The direct forcing of the second and third harmonics by the blowing velocity is relatively weak, while the excitation by nonlinear interactions of the fundamental is more substantial. This also explains their more oscillatory profiles.

Figure 10. Profiles of $Q_{n}$ in the CL for different roughness height.

Figure 11. Profiles of the leading-order chordwise velocity harmonics $\tilde {u}_{n}^{\langle 1 \rangle }$ in the CL for different roughness height.

Figure 12 $(a)$ shows the profiles of the leading-order chordwise velocity of the mean-flow distortion (2.59) in the CL. The distribution features a high shear near the critical level $y = y_{c}$ $(\tilde {Y} = 0)$ . Note that although the velocity of the distortion is $\mathit{O} (R^{-1/3})$ , the associated shear, or vorticity, is $\mathit{O} (1)$ . The mean-flow distortion is not confined because it exhibits a jump across the CL whereby driving a mean-flow distortion of the same order-of-magnitude in the main layer. The latter must be considered to determine the constant $ \tilde {u}_{0}^{\langle 1 \rangle } |_{\tilde {Y} \to \infty }$ . Figure 12 $(b)$ displays the variations with roughness height $h$ of the velocity jump (2.61) and shear jump (2.71a ) of the leading- and second-order mean-flow distortions. The magnitudes of both jumps increase with $h$ .

Figure 12. The mean-flow distortion in the CL: $(a)$ the leading-order chordwise velocity profiles for different roughness height and $(b)$ jumps of the mean-flow distortion versus roughness height.

A significant feature of the CL dynamics is that the regularised disturbance velocity, which is of $\mathit{O} (R^{-1/3})$ , significantly exceeds that in the main layer. Similar to the WL, the amplitude of the associated vorticity in the CL can be of $\mathit{O} (1)$ . We also rotate the $x$ – $z$ plane to align the $x$ -axis with the direction of $(\tilde {\alpha },\tilde {\beta })$ , where $\{ \tilde {\alpha },\tilde {\beta } \}$ , introduced in (4.1), are wavenumbers of the secondary instability mode in the distorted CL flow. The chordwise and spanwise velocity components are projected accordingly, as specified in (3.17) with $\bar {\theta }$ being replaced by $\tilde {\theta } = \arctan (\tilde {\beta } / \tilde {\alpha })$ . The vorticity along the direction of $(\tilde {\alpha },\tilde {\beta })$ can be calculated as

(3.20) \begin{equation} \tilde {\varOmega } = - \left ( \tilde {\alpha } \Big / \sqrt {\tilde {\alpha }^{2} + \tilde {\beta }^{2}} \right ) \left [ (\tilde {\beta } / \tilde {\alpha }) \left( U_{\!B}^{\prime} (y_{c}) + \tilde {U}_{\tilde {Y}} \right) - \left( W_{\!B}^{\prime} (y_{c}) + \tilde {W}_{ \tilde {Y}} \right) \right ] + {\cdots} . \end{equation}

Figure 13 displays contours of this skewed CL vorticity without the contributions from base-flow wall shears. For small $h$ , the vorticity is rather weak. However, for $h = 0.43$ (figure 13 b), a significantly distorted flow field is observed. The vorticity in the direction perpendicular to $(\tilde {\alpha },\tilde {\beta })$ controls the local secondary instability in the CL, and its contour plot for the roughness with $h = 0.43$ will be displayed later in § 5.1, where instability results are discussed.

Figure 13. Contours of the skewed CL vorticity $\tilde {\varOmega }$ without the contributions from the base-flow shear in Case I: $(a)$ $h = 0.1$ and $(b)$ $h = 0.43$ .

4.1. Secondary instability in the CL

We start with investigating spatial secondary instability in the CL. The wall-normal length scale of the local flow is $\mathit{O} (R^{-1/3})$ . The instability is taken to have characteristic wavelength of this order. Since the local velocity $(U_{\!B} (y_c), W_{\!B} (y_c) )$ is of $\mathit{O} (1)$ , the instability would have $\mathit{O} (R^{1/3})$ frequencies. With such a perception in mind, we perturb the nonlinearly distorted CL flow by a disturbance of $\mathit{O} (\epsilon _{s})$ , and express the locally perturbed flow as

(4.1) \begin{align} (u,v,w,p) &= \bigl (U_B(y_c) + R^{-1/3} U_{\!B}^{\prime}(y_c) \tilde {Y}, R^{-1} V_B(y_c), W_B(y_c) + R^{-1/3} W_{B}^{\prime}(y_c) \tilde {Y}, P_B \bigr ) \nonumber \\ & \quad + \left (R^{-1/3} \tilde {U} (x, \tilde {Y},z), R^{-2/3} \tilde {V} (x, \tilde {Y}, z), R^{-1/3} \tilde {W} (x, \tilde {Y},z), R^{-2/3} \tilde {P} (x, z) \right ) \nonumber \\ & \quad + \epsilon _{s} \left [ \left (\tilde {U}^{\langle s \rangle } (x, \tilde {Y}, z), \tilde {V}^{\langle s \rangle } (x, \tilde {Y}, z), \tilde {W}^{\langle s \rangle } (x, \tilde {Y}, z), R^{-1/3} \tilde {P}^{\langle s \rangle } (x, \tilde {Y}, z) \right )\right . \nonumber \\ & \quad \left . \times \textrm{exp} \left ( \textrm{i} R^{1/3} ( \tilde {\alpha } x + \tilde {\beta } z - \tilde {\omega } t ) \right ) + \textrm{c.c.} \right ]\!, \end{align}

where the last term represents a secondary instability mode and $\epsilon _{s} \ll 1$ . Note that the factor $R^{1/3}$ in the exponent indicates the high-frequency short-wavelength character of the instability with $\tilde {\omega }$ and $(\tilde {\alpha }, \tilde {\beta })$ representing the rescaled frequency and wavenumbers. Furthermore, we require that

(4.2) \begin{equation} \tilde {\omega } = \tilde {\alpha } U_{\!B} (y_c) + \tilde {\beta } W_{\!B} (y_c) \end{equation}

such that the secondary instability is trapped within the CL. For this spatial secondary instability, we fix the high-frequency $\tilde {\omega }$ , or equivalently the wavenumbers $\{ \tilde {\alpha }, \tilde {\beta } \}$ , to be real, and seek $\mathit{O} (1)$ modal growth in the $x$ -direction. Substituting (4.1) into the N–S (2.3), we obtain the governing equations for the spatial CL mode,

(4.3a) \begin{align} \textrm{i} \tilde {\alpha } \tilde {U}^{\langle s \rangle } + \tilde {V}_{\tilde {Y}}^{\langle s \rangle } + \textrm{i} \tilde {\beta } \tilde {W}^{\langle s \rangle } = 0, \end{align}

(4.3b) \begin{align} \tilde {\mathcal{L}} \tilde {U}^{\langle s \rangle } + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \tilde {U}^{\langle s \rangle } + (U_{\!B}^{\prime}(y_c)+\tilde {U}_{\tilde {Y}}) \tilde {V}^{\langle s \rangle } + \textrm{i} \tilde {\alpha } \tilde {P}^{\langle s \rangle } = 0, \end{align}

(4.3c) \begin{align} \tilde {\mathcal{L}} \tilde {V}^{\langle s \rangle } + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \tilde {V}^{\langle s \rangle } + \tilde {P}_{\tilde {Y}}^{\langle s \rangle } = 0, \end{align}

(4.3d) \begin{align} \tilde {\mathcal{L}} \tilde {W}^{\langle s \rangle } + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \tilde {W}^{\langle s \rangle } + (W_{B}^{\prime}(y_c)+\tilde {W}_{\tilde {Y}} ) \tilde {V}^{\langle s \rangle } + \textrm{i} \tilde {\beta } \tilde {P}^{\langle s \rangle } = 0, \end{align}

where we have introduced the operator

(4.4) \begin{equation} \tilde {\mathcal{L}} = U_B(y_c) \frac {\partial }{\partial x} + W_B(y_c) \frac {\partial }{\partial z} + \textrm{i} \bigl [ \tilde {\alpha } U_{\!B}^{\prime}(y_c) + \tilde {\beta } W_{B}^{\prime}(y_c) \bigr ] \tilde {Y}. \end{equation}

By eliminating $\tilde {U}^{\langle s \rangle }$ , $\tilde {W}^{\langle s \rangle }$ and $\tilde {P}^{\langle s \rangle }$ among the equations in (4.3), the system is reduced to a single equation for $\tilde {V}^{\langle s \rangle }$ ,

(4.5) \begin{equation} \left [ \tilde {\mathcal{L}} + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \right ] \left (\frac {\partial ^2}{\partial \tilde {Y}^2} - \tilde {k}^{2} \right ) \tilde {V}^{\langle s \rangle } - \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} )_{ \tilde {Y} \tilde {Y} } \tilde {V}^{\langle s \rangle } = 0, \end{equation}

where $\tilde {k} = \sqrt {\tilde {\alpha }^{2} + \tilde {\beta }^{2}}$ .

To formulate the eigenvalue problem for the spatial CL mode, we assume a modal dependence of $x$ and express the instability solution in the form,

(4.6) \begin{equation} ( \tilde {U}^{\langle s \rangle }, \tilde {V}^{\langle s \rangle }, \tilde {W}^{\langle s \rangle }, \tilde {P}^{\langle s \rangle } ) = ( \tilde {u}^{\langle s \rangle }, \tilde {v}^{\langle s \rangle }, \tilde {w}^{\langle s \rangle }, \tilde {p}^{\langle s \rangle } ) \textrm{exp} ( \tilde {\sigma } x ), \end{equation}

where $\tilde {\sigma } = \tilde {\sigma }_{r} + \textrm{i} \tilde {\sigma }_{i}$ is to be sought as the eigenvalue, with $\tilde {\sigma }_{r}$ representing the growth rate, which is of $\mathit{O} (1)$ . Substituting (4.6) into (4.5), we obtain the governing equation

(4.7) \begin{equation} \Bigl [ \tilde {\sigma } U_{\!B} (y_{c}) + \tilde {\mathcal{L}} + \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} ) \Bigr ] \left (\frac {\partial ^2}{\partial \tilde {Y}^2} - \tilde {k}^{2} \right ) \tilde {v}^{\langle s \rangle } - \textrm{i} ( \tilde {\alpha } \tilde {U} + \tilde {\beta } \tilde {W} )_{ \tilde {Y} \tilde {Y} } \tilde {v}^{\langle s \rangle } = 0. \end{equation}

Since the leading-order roughness-induced disturbance concentrates in the CL, (4.7) collapses to $(\partial _{\tilde {Y} \tilde {Y}}^{2} - \tilde {k}^{2} ) \tilde {v}^{\langle s \rangle } \to 0$ as $\tilde {Y} \to \pm \infty$ , which is further simplified to give the boundary conditions for the CL mode

(4.8) \begin{equation} \left ( \frac {\partial }{\partial \tilde {Y}} \mp \tilde {k} \right ) \tilde {v}^{\langle s \rangle } = 0 \quad \textrm{as} \, \, \, \tilde {Y} \to \mp \infty ; \end{equation}

these serve to exclude the non-physical unbounded (exponentially large) solutions.

The roughness-induced disturbance velocities $\tilde {U}$ and $\tilde {W}$ are expressed as Fourier series, as indicated by (2.52), and we only need to retain the leading-order terms here,

(4.9) \begin{equation} (\tilde {U}, \tilde {W}) = \sum _{n = - \infty }^{\infty } \left ( \tilde {u}_{n}^{\langle 1 \rangle } (\tilde {Y}), \tilde {w}_{n}^{\langle 1 \rangle } (\tilde {Y}) \right ) \textrm{e}^{\textrm{i} n \alpha _{w} \zeta } + {\cdots} . \end{equation}

The coefficients of the (4.7) are periodic functions of $\zeta$ , and thus according to Floquet theory we express $\tilde {v}^{\langle s \rangle } (x, \tilde {Y}, z)$ in the form of Fourier series with respect to $\zeta$ ,

(4.10) \begin{equation} \tilde {v}^{\langle s \rangle } (x, \tilde {Y}, z) = \ \sum _{n = - \infty }^{\infty } \tilde {v}_{n}^{\langle s \rangle } (\tilde {Y}) \textrm{e}^{\textrm{i} n \alpha _{w} \zeta }; \end{equation}

in general the right-hand side may have an exponential factor $\textrm{exp} (\textrm{i} q \alpha _{w} \zeta )$ with $0 \le q \le 1/2$ being a Floquet-exponent-like parameter. It turns out that the resulting extra term $\textrm{i} q \beta _{w} W_{\!B} (y_c)$ can be absorbed into $\tilde {\sigma }$ , and thus we set $q = 0$ without losing generality. Substituting (4.9) and (4.10) into (4.7) with the definition of the critical $y_{c}$ leads to the governing equation for $\tilde {v}_{n}^{\langle s \rangle }$ ,

(4.11) \begin{align} & \textrm{i} ( \tilde {\alpha } U_{\!B}^{\prime}(y_c) + \tilde {\beta } W_{B}^{\prime}(y_c)) \tilde {Y} \left (\!\frac {\textrm{d}^2}{\textrm{d} \tilde {Y}^2} - \tilde {k}^{2} \!\right ) \tilde {v}_{n}^{\langle s \rangle } + \sum _{k = -\infty }^{\infty } \textrm{i} \left [ \!\left(\tilde {\alpha } \tilde {u}_{n-k}^{\langle 1 \rangle } + \tilde {\beta } \tilde {w}_{n - k}^{\langle 1 \rangle } \!\right) \!\left (\!\frac {\textrm{d}^2}{\textrm{d} \tilde {Y}^2} - \tilde {k}^{2} \!\right ) \tilde {v}_{k}^{\langle s \rangle } \right . \nonumber \\ & \quad \left . - ( \tilde {\alpha } \tilde {u}_{n-k}^{\langle 1 \rangle } + \tilde {\beta } \tilde {w}_{n-k}^{\langle 1 \rangle } )_{ \tilde {Y} \tilde {Y} } \tilde {v}_{k}^{\langle s \rangle } \right ] = - \tilde {\sigma } U_{\!B} (y_c) \left (\frac {\textrm{d}^2}{\textrm{d} \tilde {Y}^2} - \tilde {k}^{2} \right ) \tilde {v}_{n}^{\langle s \rangle }. \end{align}

The above equations are subject to boundary conditions (4.8). Discretisation of (4.11) and (4.8) gives rise to a linear generalised eigenvalue problem,

(4.12) \begin{equation} \tilde {\boldsymbol{A}} \tilde {\boldsymbol{V}}^{\langle s \rangle } = \tilde {\sigma } \tilde {\boldsymbol{B}} \tilde {\boldsymbol{V}}^{\langle s \rangle }, \end{equation}

where the solution vector $\tilde {\boldsymbol{V}}^{\langle s \rangle }$ , the matrices $\tilde {\boldsymbol{A}}$ and $\tilde {\boldsymbol{B}}$ are defined in the supplementary material (§ S4), and further details can be found there.

For temporal secondary instability, the perturbation in (4.1) is replaced by the CL mode of the form,

(4.13) \begin{align} & \epsilon _{s} \Bigl [ \bigl (\check {U}^{\langle s \rangle } (x, \tilde {Y}, z), \check {V}^{\langle s \rangle } (x, \tilde {Y}, z), \check {W}^{\langle s \rangle } (x, \tilde {Y}, z), R^{-1/3} \check {P}^{\langle s \rangle } (x, \tilde {Y}, z) \bigr ) \nonumber \\ & \qquad \qquad \qquad \qquad \qquad \qquad \times \textrm{exp} \bigl ( \check {\sigma } t + \textrm{i} R^{1/3} ( \tilde {\alpha } x + \tilde {\beta } z - \tilde {\omega } t ) \bigr ) + \textrm{c.c.} \Bigr ], \end{align}

where the condition (4.2) is required to hold, and $\check {\sigma }$ stands as the eigenvalue, with its real part representing the temporal growth rate. By expressing $\check {V}^{\langle s \rangle } (x, \tilde {Y}, z)$ as a Fourier series with respect to $\zeta$ ,

(4.14) \begin{equation} \check {V}^{\langle s \rangle } (x, \tilde {Y}, z) = \sum _{k = - \infty }^{\infty } \check {V}_{n}^{\langle s \rangle } \textrm{e}^{\textrm{i} n \alpha _{w} \zeta }, \end{equation}

the governing equation for the Fourier component $\check {V}_{n}^{\langle s \rangle }$ can be derived and turns out to be the same as that for spatial instability, i.e. (4.11), with $\check {\sigma }$ playing exactly the role of $\tilde {\sigma } U_{\!B} (y_{c})$ , leading to the equivalence

(4.15) \begin{equation} \check {\sigma } = U_{\!B} (y_c) \tilde {\sigma }. \end{equation}

Therefore, we will not treat the temporal instability separately, and only show numerical results for spatial CL modes. As this small-scale local instability is newly identified, verification through initial-value calculations is necessary; the formulation and numerical results are presented in the supplementary material (§ S3).

4.2. Secondary instability in the WL

We now analyse instability of the roughness-distorted WL flow to small-scale perturbation. The WL has an $\mathit{O} (R^{-1/3})$ width, and we seek WL modes with characteristic wavelength of this order. Unlike the CL, the local velocity is $\mathit{O} (R^{-1/3})$ , and so the characteristic frequency of WL modes is of $\mathit{O} (1)$ .

4.2.1. Temporal modes

We first study the temporal WL mode. The instantaneous WL flow, consisting of the roughness-distorted part and the $\mathit{O} (\epsilon _{s})$ instability mode, is expressed as

(4.16) \begin{align} (u, v, w, p) &= \left(\! R^{-1/3} \bar {U}_{w} (x, \bar {Y}\!, z), R^{-2/3} \bar {V}_{w} (x, \bar {Y}\!, z), R^{-1/3} \bar {W}_{w} (x, \bar {Y}\!, z), R^{-2/3} \bar {P}_{w} (x, z)\!\right) \nonumber \\ & \quad + \epsilon _{s} \bigl [ \bigl ( \bar {U}^{\langle s \rangle } (x, \bar {Y}, z), \bar {V}^{\langle s \rangle } (x, \bar {Y}, z), \bar {W}^{\langle s \rangle } (x, \bar {Y}, z), R^{-1/3} \bar {P}^{\langle s \rangle } (x, \bar {Y}, z) \bigr ) \nonumber \\ & \quad \times \textrm{exp} \bigl ( \bar {\sigma } t + \textrm{i} R^{1/3}(\bar {\alpha } x + \bar {\beta } z ) \bigr ) + \textrm{c.c.} \bigr ], \end{align}

where the last term represents a WL mode with its scaled chordwise and spanwise wavenumbers $\bar {\alpha }$ and $\bar {\beta }$ being taken as the real parameters, and $\bar {\sigma } = \bar {\sigma }_{r} + \textrm{i} \bar {\sigma }_{i}$ is to be sought as the eigenvalue, whose real and imaginary parts represent the temporal growth rate and frequency of the WL mode, respectively. By inserting (4.16) into (2.3), the governing equations for the temporal WL mode are derived as

(4.17a) \begin{align} \textrm{i} \bar {\alpha } \bar {U}^{\langle s \rangle } + \bar {V}_{\bar {Y}}^{\langle s \rangle } + \textrm{i} \bar {\beta } \bar {W}^{\langle s \rangle } = 0, \end{align}

(4.17b) \begin{align} \bar {\sigma } \bar {U}^{\langle s \rangle } + \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right ) \bar {U}^{\langle s \rangle } + \bar {V}^{\langle s \rangle } \bar {U}_{w, \bar {Y}} + \textrm{i} \bar {\alpha } \bar {P}^{\langle s \rangle } = 0, \end{align}

(4.17c) \begin{align} \bar {\sigma } \bar {V}^{\langle s \rangle } + \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right ) \bar {V}^{\langle s \rangle } + \bar {P}_{\bar {Y}}^{\langle s \rangle } = 0, \end{align}

(4.17d) \begin{align} \bar {\sigma } \bar {W}^{\langle s \rangle } + \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right ) \bar {W}^{\langle s \rangle } + \bar {V}^{\langle s \rangle } \bar {W}_{w, \bar {Y}} + \textrm{i} \bar {\beta } \bar {P}^{\langle s \rangle } = 0. \end{align}

Equations in (4.17d ) can be reduced into a single equation for $\bar {V}^{\langle s \rangle }$ ,

(4.18) \begin{equation} \Bigl [ \bar {\sigma } + \textrm{i} ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} ) \Bigr ] \left ( \frac {\partial ^2}{\partial \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}^{\langle s \rangle } - \textrm{i} \left ( \bar {\alpha } \bar {U}_{w} + \bar {\beta } \bar {W}_{w} \right )_{\bar {Y} \bar {Y}} \bar {V}^{\langle s \rangle } = 0, \end{equation}

where $\bar {k} = \sqrt { \bar {\alpha }^{2} + \bar {\beta }^{2} }$ . Since the roughness-distorted WL flow, given by (2.25) and calculated in § 3, is expressed in the form of Fourier series of $\zeta$ , the WL mode can be expressed in the same form,

(4.19) \begin{equation} \bar {V}^{\langle s \rangle } \left (x, \bar {Y}, z \right ) = \sum _{n = -\infty }^{\infty } \bar {V}_{n}^{\langle s \rangle } \textrm{e}^{\textrm{i} n \alpha _w \zeta }. \end{equation}

Substitution of (2.25) and (4.19) into (4.18) leads to

(4.20) \begin{align} & \textrm{i} ( \bar {\alpha } \lambda _1 + \bar {\beta } \lambda _3 ) \bar {Y} \left ( \frac {\textrm{d}^2}{\textrm{d} \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}_{n}^{\langle s \rangle } \nonumber \\ & \quad + h \sum _{k = -\infty }^{\infty } \textrm{i} \biggl \{ \Bigl [ \bar {\alpha } \left ( \lambda _1 F_{n-k} + \bar {U}_{n-k} \right ) + \bar {\beta } \left ( \lambda _3 F_{n-k} + \bar {W}_{n-k} \right ) \Bigr ] \left ( \frac {\textrm{d}^2}{\textrm{d} \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}_{k}^{\langle s \rangle } \nonumber \\ & \quad - ( \bar {\alpha } \bar {U}_{n-k} + \bar {\beta } \bar {W}_{n-k} )_{\bar {Y} \bar {Y}} \bar {V}_{k}^{\langle s \rangle } \biggr \} = -\bar {\sigma } \left ( \frac {\textrm{d}^2}{\textrm{d} \bar {Y}^2} - \bar {k}^{2} \right ) \bar {V}_{n}^{\langle s \rangle }. \end{align}

The secondary instability mode is trapped within the WL, that is, it decays exponentially as $\bar {Y} \to \infty$ , as can be inferred from (4.18). However, unlike the CL mode, which attenuates exponentially beneath the CL, the WL mode has to satisfy the impermeability condition at the wall surface. Therefore, the boundary conditions for (4.20) are

(4.21a,b) \begin{align} \left ( \frac {\textrm{d}}{\textrm{d} \bar {Y}} + \bar {k} \right ) \bar {V}_{n}^{\langle s \rangle } = 0 \quad \textrm{as} \, \, \, \bar {Y} \to \infty ; \qquad \bar {V}_{n}^{\langle s \rangle } = 0 \quad \textrm{at} \, \, \, \bar {Y} = 0. \end{align}

Again, discretisation of (4.20) and (4.21) leads to a linear generalised eigenvalue problem similar to (4.12),

(4.22) \begin{equation} \bar {\boldsymbol{A}} \bar {\boldsymbol{V}}^{\langle s \rangle } = \bar {\sigma } \bar {\boldsymbol{B}}\bar {\boldsymbol{V}}^{\langle s \rangle }, \end{equation}

where the matrices and the details of the formulation are provided in the supplementary material (§ S4).

4.2.2. Spatial modes

Finally, we consider the spatial WL mode. The temporal mode in (4.16) is replaced by

(4.23) \begin{align} & \epsilon _{s} \bigl [ (\bar {\bar {U}}^{\langle s \rangle } (x, \bar {Y}, z), \bar {\bar {V}}^{\langle s \rangle } (x, \bar {Y}, z), \bar {\bar {W}}^{\langle s \rangle } (x, \bar {Y}, z), R^{- 1 / 3} \bar{\bar{P}}^{\langle s \rangle} (x, \bar {Y}, z)) \nonumber \\ & \quad \times \textrm{exp} \bigl (\textrm{i} R^{1/3} (\bar {\bar {\alpha }} x + \bar {\bar {\beta }} z ) - \textrm{i} \bar {\bar {\omega }} t \bigr ) + \textrm{c.c.} \bigr ], \end{align}

where $\bar {\bar {\beta }}$ and $\bar {\bar {\omega }}$ are taken as real parameters, and $\bar {\bar {\alpha }} = \bar {\bar {\alpha }}_{r} + \textrm{i} \bar {\bar {\alpha }}_{i}$ is the eigenvalue with $\bar {\bar {\alpha }}_{r}$ and $- \bar {\bar {\alpha }}_{i}$ representing the wavenumber and the spatial growth rate, respectively. Note that the spatial growth rate is of $\mathit{O} (R^{1/3})$ . The governing equation for $\bar {\bar {V}}^{\langle s \rangle }$ can be derived as

(4.24) \begin{equation} \Bigl [ - \bar {\bar {\omega }} + \left( \bar {\bar {\alpha }} \bar {U}_{w} + \bar {\bar {\beta }} \bar {W}_{w} \right) \Bigr ] \left ( \frac {\partial ^2}{\partial \bar {Y}^2} - \bar {\bar {k}}^{2} \right ) \bar {\bar {V}}^{\langle s \rangle } - ( \bar {\bar {\alpha }} \bar {U}_{w} + \bar {\bar {\beta }} \bar {W}_{w} )_{\bar {Y} \bar {Y}} \bar {\bar {V}}^{\langle s \rangle } = 0, \end{equation}

where $\bar {\bar {k}} = \sqrt {\bar {\bar {\alpha }}^{2} + \bar {\bar {\beta }}^{2}}$ . We express $\bar {\bar {V}}^{\langle s \rangle }$ as a Fourier series of $\zeta$ and substitute it, along with (2.25), into (4.24). The resulting system and the boundary conditions form a nonlinear eigenvalue problem, which is, after truncation and discretisation, converted into an algebraic system

(4.25) \begin{equation} \bar {\bar {\boldsymbol{M}}} \bar {\bar {\boldsymbol{V}}}^{\langle s \rangle } = 0. \end{equation}

The requirement of a non-trivial solution leads to an implicit dispersion relation for the spatial secondary instability

(4.26) \begin{equation} \mathcal{D} (\bar {\bar {\alpha }}; \bar {\bar {\beta }}, \bar {\bar {\omega }}) \equiv \textrm{det} \left ( \bar {\bar {\boldsymbol{M}}} (\bar {\bar {\alpha }}; \bar {\bar {\beta }}, \bar {\bar {\omega }}) \right ) = 0. \end{equation}

The derivation and the matrix $\bar {\bar {\boldsymbol{M}}}$ are given in the supplementary material (§ S4).

5.1. Secondary instability in the CL

For the secondary instability in the CL, we only calculate the spatial CL mode by solving the linear eigenvalue formulation (4.12) directly using an appropriate MATLAB eigenvalue solver. It is useful to note symmetry properties of the eigenvalues. First, (4.11) indicates that for each wavenumber pair $(\tilde {\alpha }, \tilde {\beta })$ , $\tilde {\sigma }$ obeys the relation

(5.1) \begin{equation} \tilde {\sigma } (-\tilde {\alpha }, -\tilde {\beta } ) = - \tilde {\sigma }(\tilde {\alpha }, \tilde {\beta }). \end{equation}

Second, if (4.7) with boundary conditions (4.8) admits an eigenvalue $\tilde {\sigma }$ with eigenfunction $\tilde {\boldsymbol{V}}^{\langle s \rangle }$ , then $- (\tilde {\sigma })_{cc}$ and $(\tilde {\boldsymbol{V}}^{\langle s \rangle })_{cc}$ also represent an eigenvalue and eigenfunction, which follows from taking complex conjugate of (4.7). For a given pair $(\tilde {\alpha }, \tilde {\beta })$ , we can calculate the most unstable eigenmode whose eigenvalue $\tilde {\sigma }$ has the largest positive real part, while $-(\tilde {\sigma })_{cc}$ is another eigenvalue corresponding to the most damped mode. On the other hand, according to (5.1), we have $-(\tilde {\sigma })_{cc} = [ \tilde {\sigma } (-\tilde {\alpha },-\tilde {\beta }) ]_{cc}$ . The latter is the most unstable eigenmode for $(-\tilde {\alpha }, -\tilde {\beta })$ . Therefore, it suffices to consider only the upper-half of the $\tilde {\alpha }$ – $\tilde {\beta }$ parameter plane, i.e. restrict the calculations to $\tilde {\beta } \gt 0$ .

Since the matrix formulation in (4.12) involves two-dimensional discretisation in the transverse direction $\tilde {Y}$ and Fourier spectral space, the resulting matrices are large, making direct eigenvalue calculation challenging. We begin by calculating the eigenvalue of (4.12) using a coarse grid with the direct method (MATLAB function eig). The results from the coarse grid are then used as initial guesses to initialise iterations in the Arnoldi method (MATLAB function eigs) to obtain more accurate eigenvalues on a refined grid.

We calculate the spatial CL mode supported by the nonlinearly distorted CLs pertinent to Cases I and II roughness, with heights $h = 0.43$ and $h = 1.135$ , respectively (which are the greatest heights for which convergent solutions for roughness-induced flows can be computed). The numerical results corresponding to Case II roughness are presented in the supplementary material (§ S2). Figure 14 presents contours of the growth rates $ \tilde {\sigma }_{r}$ in the upper-half of the $\tilde {\alpha }$ – $\tilde {\beta }$ plane. There exist two families of modes with the maximal growth rate $\tilde {\sigma }_{r} = 0.0280$ located at $(\tilde {\alpha }, \tilde {\beta }) = (-0.21,0.05)$ . It is noteworthy that the leading-order spanwise velocity of the roughness-distorted flow $\tilde {w}_{n}^{\langle 1 \rangle }$ is proportional to its chordwise counterpart $\tilde {u}_{n}^{\langle 1 \rangle }$ by the factor $- \alpha _{w} / \beta _{w}$ , as indicated by (2.54a ). This implies that when $(\tilde {\alpha }, \tilde {\beta })$ aligns with $(\alpha _{w}, \beta _{w})$ , the leading-order roughness-induced flow controlling the secondary instability is cancelled out in (4.11), which becomes degenerated. The present dominant instability no longer exists, and a weaker instability might be possible if the higher-order terms, $\tilde {u}_{n}^{\langle 2 \rangle }$ and $\tilde {w}_{n}^{\langle 2 \rangle }$ , are taken into account.

Figure 14. Contours of the growth rate $ \tilde {\sigma }_{r}$ of the spatial CL mode on the nonlinear CL flow with $h = 0.43$ . The maximal growth rate is indicated by the red cross. The red dashed lines represent the direction parallel to $(\alpha _{w}, \beta _{w})$ , along which the leading-order effect of the roughness-induced disturbance cancels out.

We now present quantitative information about the instability characteristics. In figure 15, the eigenvalues are plotted for three typical values of $\tilde {\beta }$ . For each $\tilde {\beta }$ , there exist two branches of modes, one branch with negative $\tilde {\alpha }$ featuring the dominant peak and a subdominant peak confined to small but negative $\tilde {\alpha }$ , and another branch with $\tilde {\alpha }$ mostly being positive. The two branches appear to be disconnected with the left-hand having $\tilde {\sigma }_{i} \lt 0$ , and the right-hand having $\tilde {\sigma }_{i} \gt 0$ . Furthermore, for the left-hand branch $\tilde {\sigma }_{i}$ varies almost linearly with respect to $\tilde {\alpha }$ , indicating that the mode is nearly non-dispersive. In contrast, the right-hand branch mode is dispersive. For $\tilde {\beta } = \{ 0.03, 0.05,0.07 \}$ chosen in figure 15, the roughness wavenumber-aligned $\tilde {\alpha }$ are $\{ -0.026, -0.043, -0.061 \}$ , respectively, the position of which are marked in figure 15 $(a)$ . The instability disappears when approaching these points. Initial-value calculations were performed to confirm some of eigenvalue results shown in figure 15 (see the supplementary material (§ S3)). When the spatial growth rates shown in figure 15 are converted using (4.15), the corresponding temporal growth rates are in the range of $(0, 0.025)$ , which is approximately twice that of the secondary instability of stationary cross-flow vortices at the same Reynolds number (Malik et al. Reference Malik, Li and Chang1994).

Figure 15. Instability characteristics of the spatial CL mode for different values of $\tilde {\beta }$ . The black-filled symbols on the $\tilde {\alpha }$ -axis mark the values of $\tilde {\alpha }$ aligning with the roughness wavenumber for each $\tilde {\beta }$ .

Figure 16 displays contours of the normalised eigenfunction $\tilde {\boldsymbol{V}}^{\langle s \rangle }$ of the most unstable mode for $(\tilde {\alpha }, \tilde {\beta }) = (-0.21,0.05)$ and $(\tilde {\alpha }, \tilde {\beta }) = (0.03,0.05)$ , respectively. Superposed are the contours of the $\mathit{O} (1)$ CL vorticity, defined as

(5.2) \begin{equation} \tilde {\varOmega }_{\perp } = (\tilde {\alpha } / \tilde {k}) \Bigl \{ U_{\!B}^{\prime} (y_c) + (\tilde {\beta } / \tilde {\alpha }) W_{\!B}^{\prime} (y_c) + \sum _{n = -\infty }^{\infty } \bigl [ \tilde {u}_{n}^{\langle 1 \rangle } + (\tilde {\beta } / \tilde {\alpha }) \tilde {w}_{n}^{\langle 1 \rangle } \bigr ]_{\tilde {Y}} \textrm{e}^{\textrm{i} n \alpha _{w} \zeta } \Bigr \}, \end{equation}

which represents the vorticity component projected to the direction perpendicular to $(\tilde {\alpha },\tilde {\beta })$ . The term in summation in (5.2) corresponds to the roughness-induced vorticity. It is the component that fundamentally alters the stability (despite small $\mathit{O} (R^{-1/3})$ velocity distortions), and indeed controls the secondary instability as (4.11) indicates. The CL mode with negative $\tilde {\alpha }$ appears rather local, residing primarily in a relatively small region where $\tilde {\varOmega }_{\perp }$ is negative (figure 16 a). With its short-scale carrier wave being accounted for, rapid oscillations in time and space are expected to be observed in this region. In contrast, the CL mode with positive $\tilde {\alpha }$ spreads out in the planar directions, thus exhibiting a distributed character (figure 16 b).

Figure 16. Contours of the normalised eigenfunction corresponding to the most unstable spatial CL mode (solid lines), superposed onto the eigenfunction contours are the skewed vorticity $\tilde {\varOmega }_{\perp }$ of the roughness-distorted CL flow (dashed lines and labelled with values): (a) $(\tilde {\alpha },\tilde {\beta }) = (-0.21,0.05)$ ; (b) $(\tilde {\alpha },\tilde {\beta }) = (0.03,0.05)$ .

5.2. Secondary instability in the WL

5.2.1. Temporal WL mode

For the temporal WL mode, the linear generalised eigenvalue problem (4.22) is first solved using a MATLAB eigenvalue solver. As the temporal mode obeys the same symmetries as the CL mode, we restrict our calculations to $\bar {\beta } \gt 0$ in the $\bar {\alpha }$ – $\bar {\beta }$ parameter plane.

We sought the temporal WL mode using the numerical solutions for the nonlinear roughness-distorted WL flow in Case I with $h = 0.825$ and Case II with $h = 1.135$ , which are the largest heights for which converged solutions can be computed. The numerical results for Case II are shown in the supplementary material (§ S2). Figure 17 shows contours of the growth rate $\bar {\sigma }_{r}$ in the upper-half of the $\bar {\alpha }$ – $\bar {\beta }$ plane, revealing a narrow band of modes exhibiting significant growth rates. The maximal growth rate, $\bar {\sigma }_{r} = 0.0071$ , is identified, attained at $(\bar {\alpha },\bar {\beta }) = (-0.092,0.185)$ . Figure 18 displays the instability characteristics for three different values of $\bar {\beta }$ . For each $\bar {\beta }$ , there exists a band of unstable modes, an interesting feature of which is the presence of a dominant peak and a subdominant peak; this is also shown in figure 17 $(b)$ , where there is a gap between two slanted regions. The linear variation of $\bar {\sigma }_{i}$ with $\bar {\alpha }$ indicates that the WL modes have almost the same phase speeds and are thus weakly dispersive, similar to the CL modes with $\tilde {\alpha } \lt 0$ .

Figure 17. Contours of the growth rate $\bar {\sigma }_{r}$ of the temporal WL mode (a) and local contours (b) marked by the rectangle in (a). The maximum is indicated by the red cross. A red diamond in the right panel indicates a subpeak for $\bar {\beta } = 0.18$ shown in figure 18 $(a)$ .

Figure 18. Instability characteristics of the temporal WL mode for different values of $\bar {\beta }$ .

We then plot in figure 19 contours of the normalised eigenfunction $\bar {\boldsymbol{V}}^{\langle s \rangle }$ of an unstable temporal WL mode, $(\bar {\alpha },\bar {\beta }) = (-0.10,0.20)$ . Superposed in the background are the contours of the $\mathit{O} (1)$ vorticity of the roughness-distorted WL flow, defined as

(5.3) \begin{equation} \bar {\varOmega }_{\perp } = (\bar {\alpha } / \bar {k}) \Bigl \{ \lambda _{1} + (\bar {\beta } / \bar {\alpha }) \lambda _{3} + h \sum _{n = - \infty }^{\infty } \bigl [ \bar {U}_{n} + (\bar {\beta } / \bar {\alpha }) \bar {W}_{n} \bigr ]_{\bar {Y}} \textrm{e}^{\textrm{i} n \alpha _{w} \zeta } \Bigr \}, \end{equation}

which represents the vorticity component projected to the direction perpendicular to that of $(\bar {\alpha },\bar {\beta })$ . The summation term in (5.3) corresponds to the roughness-induced vorticity in the WL, and it is this term that alters the WL stability, and indeed governs the temporal WL mode as can be seen in (4.20). Similar to the CL mode with negative $\tilde {\alpha }$ , the WL mode appears also fairly localised in its horizontal extent that is a function of the wavelength of the roughness arrays.

Figure 19. Contours of the normalised eigenfunction corresponding to an unstable temporal WL mode (solid lines), superimposed are contours of the skewed vorticity $\bar {\varOmega }_{\perp }$ of the roughness-distorted flow (dashed lines and labelled with values).

5.2.2. Spatial WL mode

The spatial WL mode amounts to a nonlinear eigenvalue problem (4.26). The eigenvalue is a root of (4.26), which needs to be obtained by an iteration process with appropriate initial guesses. The discretisation in two dimensions leads to the size of matrix $\bar {\bar {\boldsymbol{M}}}$ being too large to evaluate its determinant directly. We circumvent this difficulty by employing instead an indirect method via a rank-revealing QR factorisation of $\bar {\bar {\boldsymbol{M}}}$ (Kublanovskaya Reference Kublanovskaya1970),

(5.4) \begin{equation} \bar {\bar {\boldsymbol{M}}} \boldsymbol{P} = \boldsymbol{Q} \boldsymbol{R}, \end{equation}

with $\boldsymbol{P}$ being a permutation matrix which ensures that the diagonal entries in $\boldsymbol{R}$ form a decreasing sequence according to their magnitude so that the last diagonal entry of $\boldsymbol{R}$ , denoted as $r_{br}$ , becomes zero before any other diagonal entry does. It follows that $\textrm{det} (\bar {\bar {\boldsymbol{M}}}) = \textrm{det} (\boldsymbol{R}) = \prod \textrm{diag}(\boldsymbol{R})$ , which implies that the roots of $r_{br} = 0$ are also the roots of $\textrm{det} (\bar {\bar {\boldsymbol{M}}}) =0$ (Güttel & Tisseur Reference Güttel and Tisseur2017). Therefore, we apply Muller’s method (Muller Reference Muller1956) to

(5.5) \begin{equation} r_{br} (\bar {\bar {\alpha }}) = 0. \end{equation}

Figure 20. Comparisons between the transformed and directly calculated spatial WL mode for $\bar {\bar {\beta }} = 0.20$ .

On the other hand, a spatial instability mode can be related to a temporal mode through certain transformations, such as the well-known Gaster transformation. As the validity and accuracy of the Gaster transformation are somewhat uncertain, we apply the consistent second-order temporal–spatial transformation recently derived (Xu et al. Reference Xu, Liu, Zhang and Wu2023). We introduce $\bar {\omega } = \textrm{i} \bar {\sigma }$ to rewrite the exponent in (4.16) as

(5.6) \begin{equation} \textrm{exp} \bigl ( \bar {\sigma } t + \textrm{i} R^{1/3} (\bar {\alpha } x + \bar {\beta } z ) \bigr ) = \textrm{exp} \bigl ( \textrm{i} R^{1/3} (\bar {\alpha } x + \bar {\beta } z) - \textrm{i} \bar {\omega } t \bigr ). \end{equation}

The consistent temporal–spatial transformation of Xu et al. (Reference Xu, Liu, Zhang and Wu2023), which links the temporal and spatial growth rates for modes of the same real frequency, i.e. $\bar {\bar {\omega }} - \bar {\omega } = - \textrm{i} \bar {\omega }_{i}$ , gives the chordwise wavenumber as

(5.7) \begin{equation} \bar {\bar {\alpha }}_{0} = \bar {\alpha } + \biggl [ -1 + \biggl ( 1 - 2 \textrm{i} \bar {\omega }_{i} \frac {\textrm{d}^{2} \bar {\omega } }{\textrm{d} \bar {\alpha }^{2}} \Big / \left ( \frac {\textrm{d} \bar {\omega }}{\textrm{d} \bar {\alpha }} \right )^{2} \biggr ) ^{1/2} \biggr ] \frac {\textrm{d} \bar {\omega }}{\textrm{d} \bar {\alpha }} \Big / \left ( \frac {\textrm{d}^{2} \bar {\omega }}{\textrm{d} \bar {\alpha }^{2}} \right )\!. \end{equation}

The transformed wavenumber $\bar {\bar {\alpha }}_{0}$ may be taken as an approximation of the spatial eigenvalue, or serve as the initial guess for the iteration process, constructed using Muller’s method, to refine the eigenvalue $\bar {\bar {\alpha }}$ .

We apply the consistent transformation to temporal instability results for $\bar {\bar {\beta }} = 0.20$ . Since the temporal–spatial transformation requires information on the derivatives $\bar {\omega }$ , we need to ensure the temporal instability modes stay on a continuous branch without entangling with other branches. The transformed chordwise wavenumbers $\bar {\bar {\alpha }}_{0}$ are shown in figure 20, where comparison is made with the refined results using the iteration. The good agreement suggests that use of the transformation to temporal instability provides spatial instability characteristics of satisfactory accuracy. Spatial mode is also non-dispersive. For a typical Reynolds number $R = 10^{3}$ , the (unrescaled) spatial growth rate of WL mode, $R^{1/3} (- \bar {\alpha }_{i})$ , is approximately $0.028$ (figure 20), comparable to that of the CL mode (figure 15). We note that these are significantly greater than (approximately twice) the typical spatial growth rates of the secondary instability of stationary cross-flow vortices reported in Högberg & Henningson (Reference Högberg and Henningson1998).

Our calculations suggest that surface roughness may induce small-scale secondary instability within local regions. The present study is, to the best of our knowledge, the first to have demonstrated the existence of such small-scale secondary instabilities, although some of previous numerical calculations found modes appear to be local in their character. Specifically, for the secondary instability of stationary cross-flow vortices, Malik et al. (Reference Malik, Li, Choudhari and Chang1999) found that the $y$ -mode and $z$ -mode reside in the bulk of the boundary layer (albeit with the former being at somewhat higher positions), and both have rather high frequencies (4500 and 3500 Hz, respectively). Wassermann & Kloker (Reference Wassermann and Kloker2002) identified a high-frequency family and a low-frequency family of modes, and the former is located at appreciably higher positions than the latter, but both were deemed to be the $z$ -mode. A further computational study by Bonfigli & Kloker (Reference Bonfigli and Kloker2007) revealed two high-frequency modes, which occupy the main boundary layer and correspond to $z$ - and $y$ -modes, respectively, and both of these modes appear to be concentrated in a thin region in the wall-normal direction. A low-frequency $z$ -mode developing close to the wall was also found and it exhibits the feature of a wall mode. These modes were detected in experiments (White & Saric Reference White and Saric2005; Serpieri & Kotsonis Reference Serpieri and Kotsonis2016). Despite visual local characters, the precise local nature of the secondary instability remains uncertain. This is because in all secondary instability analyses, the primary cross-flow vortices in the entire boundary layer were taken as the base flow, and the length scale was taken to be comparable to the local boundary-layer thickness, with the boundary conditions being imposed at the wall and at the outer edge of the background boundary layer, respectively, unlike the present formulation, which specifies explicitly the short wavelength and is pertinent to the local flow field. The boundary conditions (4.8) and (4.21) of the CL and WL modes indicate that these modes are trapped within the CL and WL, respectively, unaffected by the overall boundary-layer flow.

The instability identified in this paper would be manifested as amplification of small-scale high-frequency disturbances, possibly leading to premature transition and ‘flow tripping’. In the presence of spatially extended and periodic roughness arrays as assumed, both temporal and spatial instabilities are likely to be relevant.

From the temporal instability perspective, when the roughness height exceeds a critical value, inevitable spatially extended incidental initial perturbations within the boundary layer, which in general consist of modal components, would undergo collective amplification in time, causing an otherwise laminar state turbulence. This of course does not happen when the roughness height is below the critical value for instability. The instability may thus be associated with flow tripping by roughness. We note that the phenomenon of tripping by streamwise isolated roughness has been explained in terms of instability and resulting transition of the distorted flow (Kurz & Kloker Reference Kurz and Kloker2016). Tripping by spatially distributed roughness array is another possibility, and could be explained in a similar vein.

From the perspective of spatial instability, one may envisage that the predicted small-scale high-frequency modes are excited at a certain upstream location, and their subsequent amplification could lead to onset of turbulence at a downstream position. The process of amplification and transition are likely to be similar to that associated with the secondary instability of stationary cross-flow vortices.

Given that roughness-induced distortions and the directly induced small-scale instability resemble primary cross-flow vortices and their secondary instability, respectively, we think that experimental investigations of the former would, in principle, be similar to those of the latter. In addition to constructing appropriate patterned roughness arrays, a possible challenge is to avoid or minimise primary cross-flow vortices. This might be achieved by deploying arrays of sufficient height in a region relatively close to the leading edge so that the roughness-induced distortions and the directly induced small-scale instability dominate.

Unfortunately, to the best of our knowledge, there has not been any experimental or direct numerical investigation aimed specifically at the local (CL) structures of the distorted boundary-layer flows by periodically distributed roughness arrays, or the resulting small-scale instability. Nevertheless, we note that there have been works on using regularly distributed roughness to suppress primary cross-flow vortices and their secondary instability (Wassermann & Kloker Reference Wassermann and Kloker2002; Saric et al. Reference Saric, West, Tufts and Reed2019; Ide, Hirota & Tokugawa Reference Ide, Hirota and Tokugawa2021; Zoppini et al. Reference Zoppini, Michelis, Ragni and Kotsonis2022; Suzuki et al. Reference Suzuki, Yakeno, Konishi, Tokugawa, Hirota, Takami and Obayashi2024). Those experimental and computational techniques could be adapted to investigate/check the scenarios/findings presented here. On the other hand, we anticipate that our mathematical framework as well as theoretical results would prompt and guide future experimental and numerical investigations into this mechanism.

Of close relevance to our work are the recent studies by Hirota and his colleagues. Ide et al. (Reference Ide, Hirota and Tokugawa2021) investigated the transition-delay effect of sinusoidal roughness elements (SREs) using DNS. Unlike discrete roughness elements, these roughness arrays were distributed sinusoidally along a direction inclined at an angle to the leading edge, in which sense rather similar to the form considered in the present paper, but their overall height is modulated by a chordwise Gaussian envelope with its maximum at a distance downstream. Intriguingly, the inclination is chosen to nearly overlap the streamline at the outer edge of the boundary layer. The DNS results suggest that compared with discrete roughness elements, SREs are more viable and effective because with greater heights they do not cause tripping (i.e. direct transition to turbulence) while exciting in their wake the desired cross-flow vortices of appropriate amplitude (the so-called control mode) to inhibit the targeted most dangerous cross-flow vortices and delay the onset of the secondary instability. The effectiveness and advantage of SREs were confirmed by experiments (Suzuki et al. Reference Suzuki, Yakeno, Konishi, Tokugawa, Hirota, Takami and Obayashi2024) while detailed mechanisms remain elusive. Hirota, Ide & Hattori (Reference Hirota, Ide and Hattori2024) further investigated the role of the inclination direction of roughness arrays, showing that the desired direction must be close to that of the inviscid free stream streamline, beyond which secondary instability takes place in the immediate wake of the roughness arrays causing rapid transition. From the mathematical viewpoint, this direction is special for (at least) two reasons. Firstly, it pertains to a distinguished cross-flow instability regime, in which stationary vortices are ‘free-streamline aligned’ and of viscous wall-mode type (Choudhari Reference Choudhari1995). Secondly, the boundary-layer response to roughness arrays with this orientation does not have a CL in the bulk. As a result, the local CL instability as described in our paper does not arise, and one might thus speculate that this may be the reason why such SREs did not cause tripping. The present theory does not apply to this case, and a new formulation is required to characterise the roughness-induced distortion and possible secondary instability. Furthermore, mathematical theories could be developed for chordwise modulated arrays employed in experiments so that the excitation of the control mode, and the interactions among the targeted incoming most dangerous mode, as well as the roughness-induced distortion could be investigated.

We believe that such local instabilities are rather common and play a pivotal role in boundary-layer transition. In particular, we believe that secondary instabilities leading to breakdown of primary stationary and travelling cross-flow vortices (eigenmodes) are likely to be of local nature, and the present asymptotic approach may be adapted to provide a firm answer as well. Moreover, in fully developed shear turbulence, such small-scale instabilities that are local in space may operate to facilitate direct energy transfer from large-scale structures to much smaller-scale fluctuations, which is non-local in spectral space. Such a process may supplement the continuous energy cascade across scales.