# Stability of the solitary wave boundary layer subject to finite amplitude
disturbances

Asim Önder 1<EMAIL_ADDRESS>Philip L.-F. Liu 1,2,3  1 Department of
Civil and Environmental Engineering, National University of Singapore,
Singapore 117576, Singapore 2 School of Civil and Environmental Engineering,
Cornell University, Ithaca, NY 14850, USA 3 Institute of Hydrological and
Oceanic Sciences, National Central University, Jhongli, Taoyuan, 320, Taiwan

###### Abstract

The stability and transition in the bottom boundary layer under a solitary
wave are analysed in the presence of finite amplitude disturbances. First, the
receptivity of the boundary layer is investigated using a linear input-output
analysis, in which the environment noise is modelled as distributed body
forces. The most dangerous perturbations in a time frame until flow reversal
are found to be arranged as counter-rotating streamwise-constant rollers. One
of these roller configurations is then selected and deployed to nonlinear
equations, and streaks of various amplitudes are generated via lift-up
mechanism. By means of secondary stability analysis and direct numerical
simulations, the dual role of streaks in the boundary-layer transition is
shown. When the amplitude of streaks remains moderate, these elongated
features remain stable until the adverse-pressure-gradient stage and have a
dampening effect on the instabilities developing thereafter. In contrast, when
the low-speed streaks reach high amplitudes exceeding 15% of free-stream
velocity at the respective phase, they become highly unstable to secondary
sinuous modes in the outer shear layers. Consequently, a subcritical
transition to turbulence, i.e., bypass transition, can be already initiated in
the favourable-pressure-gradient region ahead of the wave crest.

##  1 Introduction

Solitary waves are long waves of permanent form, which induce approximately
constant velocity in the water column (Munk, 1949). They are subject to
friction in the thin boundary layers developing at the free surface and at the
sea bottom. The free-surface boundary layer is usually weak (Klettner & Eames,
2012), and is negligible. On the other hand, the bottom boundary layer is of
prominent importance, as it hosts the hydrodynamic processes driving the
sediment motion and energy dissipation. In the most basic setting of wave
propagating in a constant depth over a smooth bottom, the bottom boundary
layer consists of regions of favourable and adverse pressure gradient (FPG and
APG) located ahead and behind of the wave crest respectively, cf. figure 1.
The boundary layer flow has tendency to remain laminar in the FPG region.
Behind the wave crest, the APG gives rise to an inflectional velocity profile
(Liu et al., 2007), cf. velocity profiles in figure 1a, and the boundary layer
becomes linearly unstable (Blondeaux et al., 2012; Verschaeve & Pedersen,
2014; Sadek et al., 2015). Experimental (Sumer et al., 2010) and numerical
(Vittori & Blondeaux, 2008; Ozdemir et al., 2013) models of solitary-wave
boundary layer (SWBL) have shown that the inflectional instability leads to
regularly spaced spanwise-oriented vortex rollers, which can break down to
small-scale turbulence in higher wave amplitudes.

Figure 1:  Conceptual sketches showing two main paths of transition to
turbulence in the bottom boundary layer under a solitary wave. Scales in the
boundary layer are exaggerated for clarity. Laminar velocity profiles are
plotted until the onset of transition. (a) Orderly route to transition via
two-dimensional modal instabilities initiated by the inflectional velocity
profile. (b) Bypass transition initiated by the receptivity of boundary layer
to finite amplitude ambient disturbances (dashed-dotted curls). The
instability is three-dimensional and of stochastic nature.

Linearly stable base flow in the FPG region does not preclude the onset of
transition in this region. Finite amplitude external perturbations such as
breaking-wave turbulence, sound or small-scale bedforms can lead to
significant growth in finite times and yield secondary base states that can be
unstable. Such subcritical transition can take place in the FPG region before
the arrival of the wave, and is analogous to bypass transition in zero-
pressure gradient (ZPG) boundary layers (Morkovin, 1969), as the modal
instability is bypassed by another noise-induced mechanism. This alternative
transition scenario is depicted for SWBL in figure 1b. Unlike the orderly
transition, whose initiation is often described simply by a critical Reynolds
number, bypass transition is a complicated problem depending on the amplitude,
frequency and type of external perturbations in addition to Reynolds number. A
recent review on the phenomenology of bypass transition can be found in Durbin
(2017). All bypass scenarios require a receptivity stage, where external
perturbations are modified and amplified by the boundary layer, and a
breakdown stage, where the most amplified modes become unstable and break into
turbulence. The receptivity stage is dominated by streamwise-oriented vortices
and streaks. The former are the surviving modes from rapid shear distortion
(Phillips, 1969) and the latter are elongated streamwise-momentum modes
produced by the former via stirring the base flow, a process known as lift-up
mechanism (Landahl, 1980; Brandt, 2014). The initial stages of streak
amplification can be linked mathematically to the nonnormality of the
linarized Navier-Stokes operator (Butler & Farrell, 1992; Trefethen et al.,
1993).

The breakdown stage is characterized by secondary instabilities and resultant
turbulent spots. Two scenarios have been observed depending on the amplitude
of streaks. When the environment forcing is strong, the lift-up mechanism can
generate highly elevated low-speed streaks acting like strong wake
perturbations on their environment. These protruding layers are susceptible to
wake-like instabilities driven by spanwise shear (Waleffe, 1995; Andersson et
al., 2001; Vaughan & Zaki, 2011). Consequently, they develop rapidly growing
sinuous undulations, and break down into turbulent spots (Matsubara &
Alfredsson, 2001; Jacobs & Durbin, 2001; Hernon et al., 2007). In contrast,
when the environment forcing is modest, the streaks are weaker and remain
confined to the near-wall region. In this case, instabilities can occur on
vertical shear layers that are slightly modulated by streaks. These
instabilities are observed to have reduced growth rates compared to reference
instabilities (Tollmien–Schlichting (TS) waves), cf. Cossu & Brandt (2004) and
Liu et al. (2008b). Therefore, introducing moderate-amplitude streaks to the
boundary layer can delay transition point to turbulence (Cossu & Brandt,
2002). Vaughan & Zaki (2011) named the two streak instabilities after the
location of their respective critical layers and called them “outer” and
“inner” modes. We will adapt a similar terminology for the streak
instabilities in the present work. The final step of breakdown stage follows
the same path for both inner and outer modes, i.e., turbulent spots at
different locations grow and amalgamate (e.g. Narasimha (1985)), and finally,
turbulent boundary layer sets in.

Unlike flat-plate boundary layers, experimental and numerical evidence for
bypass transition in SWBLs are sparse. Using direct numerical simulations
(DNS), Ozdemir et al. (2013) examined the effect of the perturbation amplitude
by seeding random noise of varying magnitudes (between $1-20\%$ of maximum
free-stream velocity) to initial conditions before the arrival of the solitary
wave. The cases with $5\%$ or more noise and $Re_{\delta}\geq 1500$, where
Reynolds number is defined using Stokes length and maximum free-stream
velocity (cf. § 2 for details), showed an initial energy amplification inside
the boundary layer lasting until another more rapidly growing amplification
mechanism takes over in the APG stage after flow reversal. They speculated
that this early perturbation growth should be due to a nonlinear viscous
instability, as it takes place in the FPG stage where velocity profiles do not
contain any inflection point, a necessary condition for inviscid instability.
However, it is more likely that the amplification is due to linear transient
nonnormal growth. Indeed, Verschaeve et al. (2017) found a strong linear
nonnormal growth in the FPG stage of SWBL if the initial perturbations are
organized as streamwise rollers. These rollers then amplify streaks by the
lift-up mechanism with a maximum growth proportional to the square of Reynolds
number. Later in the APG stage, the nonnormal growth of streaks are overtaken
by the nonnormal growth of two-dimensional spanwise modes, as the maximum
growth of these modes has an exponential scaling in Reynolds number. The
analysis of Verschaeve et al. (2017) provided conceptual insights for the
overtaking of another growth mechanism in the APG stage in Ozdemir et al.
(2013). However, there is some discrepancy between DNS and transient-growth
analysis at which critical Reynolds number the nonnormal two-dimensional modes
begins to exert dominance. This is possibly related to nonlinear effects,
which are not considered in the study of Verschaeve et al. (2017). The
transition process in Ozdemir et al. (2013) was initiated by regularly spaced
vortex rollers in all cases. Although bypass transition via streak breakdown
was not observed, the secondary mechanisms in transition were sensitive to the
level of initially seeded perturbation suggesting a sensitivity to the
presence of the streaks. Low-noise cases followed a transition path
reminiscent of free-shear layers. In contrast, in high-noise cases, where the
streaks should be strongest, vortex rollers broke into $\Lambda$-shaped
vortices, hence the transition was reminiscent of K-type transition in flat-
plate boundary layers (Klebanoff et al., 1962).

Sumer et al. (2010) simulated a SWBL in an oscillatory water tunnel and
observed turbulent spots in a flow regime starting at $\Rey_{\delta}=1000$.
These sporadic features spread to earlier phases with increasing Reynolds
number. At the highest Reynolds number achieved ($\Rey_{\delta}=2000$), they
were already observed at the end of FPG stage. Turbulent spots grew and merged
but did not lead to a complete breakdown in any of the cases. They coexisted
concurrently with vortex rollers, which emerged later in laminar areas
surrounding turbulent spots, cf. video 3 in supplemental materials of Sumer et
al. (2010). The precursor structures to these turbulent spots are not clear.
It is difficult to discern streamwise streaks in the visualizations and
supplemental movies of Sumer et al. (2010). However, in a prequel paper
(Carstensen et al., 2010), where an oscillatory flow was investigated in the
same apparatus, it was clearly shown that streamwise streaks are precursors to
turbulent spots in transitional oscillatory boundary layers. It is possible
that the same facility noise induced streaks in both periodic and solitary
wave boundary layers. To date, the streak breakdown in SWBLs is explicitly
shown only in the DNS study of Sadek (2015), where a bypass scenario is
initiated by seeding optimal streamwise-constant perturbations and some
localized secondary perturbations towards the end of the FPG stage. The
injected streamwise streaks became unstable and the breakdown to small-scale
turbulence took place in the APG stage.

Studies by Ozdemir et al. (2013) and Verschaeve et al. (2017) imply that the
FPG stage of SWBL is receptive to environment perturbations and can respond by
developing streaks. There is some experimental evidence that these streaks
might breakdown into turbulent spots (Sumer et al., 2010) or modify the
secondary modes of transition when they have modest energy (Ozdemir et al.,
2013). Furthermore, we anticipate that modest-amplitude streaks can have a
stabilizing effect on the instabilities developing in the APG stage as in
flat-plate boundary layers (Cossu & Brandt, 2004). There is a need for a
systematic study to determine the quantitative and qualitative extent of these
effects. In particular, the receptivity and breakdown stages of bypass
transition in SWBLs have to be characterised in more detail. The present study
is an effort in this direction.

The disturbances in the sea environment continuously force the wave boundary
layer in a random fashion. To this end, the flow structures dominating the
early landscape in the boundary layer are induced by the external
perturbations to which the boundary layer shows the strongest response. These
perturbations can be identified in an optimisation framework using a system
perspective, where external disturbances provide the input and the boundary-
layer response corresponds to the output. A convenient approach to model the
external perturbations is using body forces of stochastic or deterministic
nature. In this context, a deterministic perturbation model allows a more
controllable approach, where the frequency and spatial distribution of body
forces can be specified. Assuming perturbations have small amplitude, a linear
approach can be utilized and response to all possible disturbances in the form
of Fourier modes can be investigated. In a pioneering work using this model,
Jovanović & Bamieh (2005) studied the linear response of the plane Pouseille
flow, and identified counter-rotating streamwise-constant rollers as the most
“dangerous” perturbation delivering the largest amplification per energy
input. The rollers induced energetic streamwise-constant streaks via the lift-
up mechanism. This study showed that the linear input-output framework can
capture the essence of receptivity stage despite its simplicity. We will use a
similar approach as a starting point in the present analysis and study the
receptivity of SWBL. We obtain input-output configurations in the form of
streak-roller systems similar to Jovanović & Bamieh (2005). Subsequently, the
breakdown stage is investigated with a combination of linear secondary
stability analysis and fully nonlinear numerical simulations triggered with
finite amplitude perturbations. We quantify the critical perturbation levels
leading to breakdown of streaks via inner and outer secondary instabilities.
Interesting results are found implying a dual role for the streaks. Weak to
moderate amplitude streaks and associated inner instabilities are shown to
have a stabilizing effect on the boundary layer, whereas higher amplitude
streaks can lead to an early bifurcation to an unstable branch already in the
FPG stage.

The manuscript is organized as follows. In §2, we briefly introduce the SWBL
in a temporal setting. Subsequently, the linear input-output framework is
described in §3, where the linear flow response is also discussed. We then
select a suitable excitation configuration and embed it to nonlinear governing
equations in §4 to obtain streaky SWBLs, which are the flow response to finite
amplitude excitation. The nonlinear flow response represents the secondary
flow states that are amenable to linear instabilities. The character and phase
of these instabilities are analysed in §5 using a linear secondary stability
analysis based on a quasi-static assumption. In §6, the growth and breakdown
of streaks is studied using nonlinear DNS. The objective of this section is
the validation of quasi-static assumption and the determination of breakdown
thresholds. Finally, in §7, the results are summed up, and implications of the
analysis are discussed with some outlook for future work.

##  2 Problem formulation

We consider a small-amplitude solitary wave with a wave height $H^{*}$
propagating over a constant depth $h^{*}$. In this work, the physical
quantities with asterisk are dimensional quantities. The problem is defined in
a Cartesian coordinate system, where $x^{*}$ is the direction of wave
propagation (also called streamwise direction), $y^{*}$ is the spanwise
direction parallel to wave crest, and $z^{*}$ is the vertical direction
extending from the bed upwards. The velocity components associated with these
directions are $u^{*}$, $v^{*}$, and $w^{*}$, respectively. The leading order
solution of solitary wave profile is given by (Grimshaw, 1970)

$\frac{\eta_{w}^{*}}{h^{*}}=1+\frac{H^{*}}{h^{*}}\mathrm{sech}^{2}\left(\sqrt{\frac{3H^{*}}{4h^{*3}}}(x^{*}-c_{w}^{*}t^{*})\right),$
(1)

where $c_{w}^{*}=\sqrt{g^{*}h^{*}}$ is the wave speed with $g^{*}$ being the
gravitational acceleration. Furthermore, the irrotational streamwise velocity,
which is constant in the water column, is given by

$u_{0}^{*}(x^{*},t^{*})=U_{0m}^{*}\mathrm{sech}^{2}\left(\sqrt{\frac{3H^{*}}{4h^{*3}}}(x^{*}-c_{w}^{*}t^{*})\right),$
(2)

with the maximum velocity

$U_{0m}^{*}=\sqrt{g^{*}h^{*}}\frac{H^{*}}{h^{*}}.$ (3)

Adjacent to the bed, there is a bottom boundary layer, in which the streamwise
velocity $u^{*}$ has also a rotational velocity component $u_{r}^{*}$. The
thickness of the boundary layer is usually much smaller than the water depth,
cf. Appendix A in Sumer et al. (2010). Therefore, streamwise variations in the
boundary layer can be considered negligible compared to temporal and vertical
variations. These assumptions allow a local parallel formulation of the bottom
boundary layer, where the flow is assumed homogeneous in horizontal
directions. At a fixed point, the irrotational velocity at the bottom (free-
stream velocity hereafter) depends now only on time and reads as follows

$u_{0}^{*}(t^{*})=U_{0m}^{*}\mathrm{sech}^{2}(-\omega_{w}^{*}t^{*}),$ (4)

where

$\omega_{w}^{*}=\sqrt{\frac{3g^{*}H^{*}}{4h^{*2}}}$ (5)

is the effective wave frequency. Using the wave frequency and kinematic
viscosity the Stokes length is defined

$\delta_{s}^{*}=\sqrt{{2\nu^{*}}/{\omega_{w}^{*}}}$ (6)

as the boundary-layer scale of the problem. Equation (4) neglects the first
and higher-order terms in ${H^{*}}/{h^{*}}$. Therefore, the model is relevant
only for ${H^{*}}/{h^{*}}\rightarrow 0$. Vittori & Blondeaux (2011) employed a
less restrictive model, in which first- and second-order terms in
${H^{*}}/{h^{*}}$ are also included. The reader is referred to their work for
the effect of the wave height on the boundary layer transition under a
solitary wave. Assuming ${H^{*}}/{h^{*}}\rightarrow 0$ and
${\delta^{*}}/{h^{*}}\rightarrow 0$ provides the advantage of reducing the
parameter space of the problem to one, the Reynolds number:

$\Rey_{\delta}=\frac{U_{0m}^{*}\delta_{s}^{*}}{\nu^{*}}.$ (7)

The Stokes length is now the only relevant length scale of the problem.

We introduce the following dimensionless velocity fields, spatial coordinates,
time and pressure, respectively,

$\mbox{\boldmath$u$}=\mbox{\boldmath$u$}^{*}/U_{0m}^{*};~{}~{}\mbox{\boldmath$x$}=\mbox{\boldmath$x$}^{*}/\delta_{s};~{}~{}t=t^{*}\omega_{w}^{*};~{}~{}p=p^{*}/\rho^{*}U_{0m}^{*2}.$
(8)

The momentum equation for the irrotational streamwise velocity at the bottom
can be expressed in local temporal frame approximation as

$\frac{2}{\Rey_{\delta}}\frac{\mathrm{d}u_{0}}{\mathrm{d}t}=-\frac{\partial
p_{0}}{\partial x},$ (9)

The free-stream pressure gradient is calculated using (4) and (9) and reads

$-\frac{\partial p_{0}}{\partial
x}=\frac{4}{\Rey_{\delta}}\mathrm{sech}^{2}(-t)\mathrm{tanh}(-t).$ (10)

The non-dimensional pressure gradient and free-stream velocity are plotted in
figures 2a and 2b, respectively. The overall balance of streamwise momentum in
the laminar SWBL is given by

$\left(\frac{\partial}{\partial t}-\frac{1}{2}\frac{\partial^{2}}{\partial
z^{2}}\right)U=2\mathrm{sech}^{2}(-t)\mathrm{tanh}(-t),$ (11)

where $U=u_{r}+u_{0}$ is the total laminar velocity containing both rotational
and irrotational components. Equation (11) is supplemented with the boundary
conditions $U(z=0,t)=0$ and $U(z\rightarrow\infty,t)=u_{0}$, and we specifiy
the initial condition $U(z,-\infty)=0$. The solution of (11) is shown in
figure 2c. This approximate model of SWBL is theoretically solved (Liu &
Orfila, 2004) and adapted in experimental (Sumer et al., 2010; Tanaka et al.,
2012) and numerical (Ozdemir et al., 2013; Sadek et al., 2015; Verschaeve et
al., 2017) studies.

Figure 2: (a) Free-stream pressure gradient (10); (b) Free-stream velocity
(4); (c) Vertical profiles of streamwise velocity in a laminar temporal
solitary-wave boundary layer. The time axis is scaled with $\pi/6$.

##  3 Optimal disturbances and the flow response

The time-dependent streamwise velocity $U(z,t)$ in (11) is the base state of
the problem, which is continuously forced by external perturbations present in
the environment. A convenient approach to study the effect of these
perturbations is to model them as body forces. In the present parallel flow
model, the forcing fields can be defined as the sum of Fourier components

$[f_{u},f_{v},f_{w}](x,y,z,t)=\iiint_{-\infty}^{\infty}[\hat{f}_{u},\hat{f}_{v},\hat{f}_{w}](z)\mathrm{e}^{\mathrm{i}(\alpha
x+\beta y+\omega_{f}t)}\mathrm{d}\alpha\mathrm{d}\beta d\omega_{f},$ (12)

where $\alpha$ and $\beta$ are streamwise and spanwise wavenumbers,
respectively, and $\omega_{f}$ is the frequency. If we consider small-
amplitude perturbations, then the flow response to each Fourier component can
be studied independently. In this linear regime, the most dangerous flow
scenarios can be initiated by finding the Fourier modes inducing the strongest
flow response. The objective of this section is to find these Fourier modes
using an optimization framework and to analyse the corresponding flow
response. In § 3.1, we introduce the optimization problem and the adjoint
method to find its solution. Subsequently, the optimal input-output
configurations and their scalings are discussed in § 3.2.

###  3.1 Methodology

We apply a Fourier ansatz in the homogeneous $x$ and $y$ directions for the
perturbation velocity and pressure

$[\tilde{u},\tilde{v},\tilde{w},\tilde{p}](x,y,z,t)=\mathrm{Re}\\{[\hat{u},\hat{v},\hat{w},\hat{p}](z,t)\mathrm{e}^{\mathrm{i}(\alpha
x+\beta y)}\\}.$ (13)

In the linear regime, each Fourier mode is excited by the corresponding
harmonic force at the same spatial wavenumber

$[f_{u},f_{v},f_{w}](x,y,z,t)=\mathrm{Re}\\{[\hat{f}_{u},\hat{f}_{v},\hat{f}_{w}](z)\mathrm{e}^{\mathrm{i}(\alpha
x+\beta y+\omega_{f}t)}\\}.$ (14)

In the present parallel flow model, the perturbation dynamics can be
conveniently studied using the forced versions of the Orr-Sommerfeld and
Squire (OSS) equations (Schmid & Brandt, 2014). To this end, two different
excitation regimes can be defined. First, when $t\ll-\pi$, the base flow is
vanishingly small, and the forcing brings the stagnant flow to a periodic
state, which is given by the set

$\displaystyle\frac{1}{\Rey_{\delta}}\left[\left(2\mathrm{i}\omega_{f}-\hat{\Delta}\right)\hat{\Delta}\right]\hat{w}_{o}(z)$
$\displaystyle=\hat{g}_{w}(z),$ (15)
$\displaystyle\frac{1}{\Rey_{\delta}}\left(2\mathrm{i}\omega_{f}-\hat{\Delta}\right)\hat{\eta}_{o}(z)$
$\displaystyle=\hat{g}_{\eta}(z),$ (16)
$\displaystyle\hat{w}_{o}(0)=\frac{\partial\hat{w}_{o}}{\partial
z}(0)=\hat{\eta}_{o}(0)$ $\displaystyle=0,$ (17)
$\displaystyle\hat{w}_{o}(z\rightarrow\infty)=\frac{\partial\hat{w}_{o}}{\partial
z}(z\rightarrow\infty)=\hat{\eta}_{o}(z\rightarrow\infty)$ $\displaystyle=0,$
(18)

where $\hat{w}$ and $\hat{\eta}$ are the vertical velocity and vertical
vorticity modes associated with the wavenumber pair ($\alpha,\beta$),
$\hat{\Delta}=k^{2}-\partial^{2}/\partial z^{2}$ represents the semi-
discretised Laplacian operator, $k^{2}=\alpha^{2}+\beta^{2}$, and
$\hat{g}_{w}$ and $\hat{g}_{\eta}$ are the external driving forces containing
the control variables
$\mbox{\boldmath$\hat{f}$}=(\hat{f}_{u},\hat{f}_{v},\hat{f}_{w})$

$\displaystyle\hat{g}_{w}$
$\displaystyle=-\mathrm{i}\alpha\frac{\partial\hat{f}_{u}}{\partial
z}-\mathrm{i}\beta\frac{\partial\hat{f}_{v}}{\partial z}-k^{2}\hat{f}_{w},$
(19) $\displaystyle\hat{g}_{\eta}$
$\displaystyle=\mathrm{i}\beta\hat{f}_{u}-\mathrm{i}\alpha\hat{f}_{v}.$ (20)

When the wave arrives, the flow is no longer periodic, and the perturbation
equations during the wave event becomes

$\displaystyle\left[\left(\frac{2}{\Rey_{\delta}}\frac{\partial}{\partial
t}+\mathrm{i}\alpha
U-\frac{1}{\Rey_{\delta}}\hat{\Delta}\right)\hat{\Delta}-\mathrm{i}\alpha\frac{\partial^{2}{U}}{\partial
z^{2}}\right]\hat{w}(z,t)$
$\displaystyle=\hat{g}_{w}(z)\mathrm{e}^{\mathrm{i}\omega_{f}t},$ (21)
$\displaystyle\left(\frac{2}{\Rey_{\delta}}\frac{\partial}{\partial
t}+\mathrm{i}\alpha
U-\frac{1}{\Rey_{\delta}}\hat{\Delta}\right)\hat{\eta}(z,t)+\mathrm{i}\beta\frac{\partial
U}{\partial z}\hat{w}(z,t)$
$\displaystyle=\hat{g}_{\eta}(z)\mathrm{e}^{\mathrm{i}\omega_{f}t},$ (22)
$\displaystyle\hat{w}(z,-\infty)=\hat{w}_{o}(z);~{}~{}~{}~{}\hat{\eta}(z,-\infty)$
$\displaystyle=\hat{\eta}_{o}(z),$ (23)
$\displaystyle\hat{w}(0,t)=\frac{\partial\hat{w}}{\partial
z}(0,t)=\hat{\eta}(0,t)$ $\displaystyle=0,$ (24)
$\displaystyle\hat{w}(z\rightarrow\infty,t)=\frac{\partial\hat{w}}{\partial
z}(z\rightarrow\infty,t)=\hat{\eta}(z\rightarrow\infty,t)$ $\displaystyle=0.$
(25)

The following compact notation is used for periodic and temporal OSS system of
equations

$L_{o}\mbox{\boldmath$\hat{q}_{o}$}=C\mbox{\boldmath$\hat{f}$},~{}~{}~{}~{}L(t)\mbox{\boldmath$\hat{q}$}=C\mbox{\boldmath$\hat{f}$}\mathrm{e}^{\mathrm{i}\omega_{f}t}$
(26)

where $\hat{\mbox{\boldmath$q$}}_{o}=[\hat{w}_{o},\hat{\eta}_{o}]$ and
$\hat{\mbox{\boldmath$q$}}=[\hat{w},\hat{\eta}]$.

The response of the flow to an excitation at a wavenumber pair
$(\alpha,\beta)$ is measured by the perturbation kinetic energy

$E(\mbox{\boldmath$\hat{u}$}):=\frac{1}{2}\int_{0}^{\infty}(|\hat{u}|^{2}+|\hat{v}|^{2}+|\hat{w}|^{2})\mathrm{d}z.$
(27)

We can express $\hat{u}$ and $\hat{v}$ in terms of $\hat{w}$ and $\hat{\eta}$,
and obtain (Schmid & Henningson, 2001)

$\displaystyle
E(\hat{\mbox{\boldmath$q$}}):=\frac{1}{2k^{2}}\int_{0}^{\infty}(k^{2}|\hat{w}|^{2}+|\frac{\partial\hat{w}}{\partial
z}|^{2}+|\hat{\eta}|^{2})\mathrm{d}z.$ (28)

We look for the most dangerous perturbations initiating the strongest response
in the linear SWBL. This is equivalent to finding an optimal control
$\mbox{\boldmath$\hat{f}$}^{opt}(z;\alpha,\beta,\omega_{f},T_{f},\Rey_{\delta})$,
which yields the maximum energy amplification at a terminal time $t=T_{f}$ per
initial energy input $E(\mbox{\boldmath$\hat{q}_{o}$})$. This is found by
solving the constrained optimization problem

$\displaystyle\begin{split}G_{f}(\alpha,\beta,\omega_{f},T_{f},\Rey_{\delta})&:=\max\limits_{\mbox{\boldmath$\hat{f}$}}\frac{E(\mbox{\boldmath$\hat{q}$}(T_{f}))}{E(\mbox{\boldmath$\hat{q}_{o}$})}\\\
&s.t.\\\ L_{o}\mbox{\boldmath$\hat{q}$}_{o}&=C\mbox{\boldmath$\hat{f}$},\\\
L(t)\mbox{\boldmath$\hat{q}$}&=C\mbox{\boldmath$\hat{f}$}\mathrm{e}^{\mathrm{i}\omega_{f}t},\\\
\|\mbox{\boldmath$\hat{f}$}\|&=1\end{split}$ (29)

where $G_{f}$ is the largest response or gain. The optimization problem (29)
is subject to constraints in the form of periodic and transient OSS systems,
and to an additional normalization constraint, which ensures the forcing
energy is unity. This optimal control analysis is closely related to the
optimal transient-growth analysis of Verschaeve et al. (2017) but differs in
control variables, i.e., instead of the growth of initial perturbations, the
response to external forcing is measured.

The optimization problem (29) is solved using an adjoint approach (Luchini &
Bottaro, 2014). In this method, a Lagrangian functional is assigned to the
optimization problem, and the optimality conditions are derived from the
stationary point of the Lagrangian, cf. Appendix A for details. To this end,
the gain $G_{f}$ is maximum when the flow is forced by the optimal forcing
configuration $\mbox{\boldmath$\hat{f}$}^{opt}$ satisfying

$\displaystyle L_{o}\mbox{\boldmath$\hat{q}$}_{o}$
$\displaystyle=C\mbox{\boldmath$\hat{f}$}^{opt},$ $\displaystyle
L(t)\mbox{\boldmath$\hat{q}$}$
$\displaystyle=C\mbox{\boldmath$\hat{f}$}^{opt}\mathrm{e}^{\mathrm{i}\omega_{f}t},$
$\displaystyle L^{+}(t)\mbox{\boldmath$\hat{q}$}^{+}$ $\displaystyle=0,$ (30)
$\displaystyle\|\mbox{\boldmath$\hat{f}$}^{opt}\|$ $\displaystyle=1,$
$\displaystyle\mbox{\boldmath$\hat{f}$}^{opt}$
$\displaystyle=\rho(\mbox{\boldmath$\hat{q}$}^{+})$ (31)

Equation (30) represents the following adjoint Orr–Sommerfeld and Squire
equations

$\displaystyle\left[\left(\frac{2}{\Rey_{\delta}}\frac{\partial}{\partial
t}+\mathrm{i}\alpha
U-\frac{1}{\Rey_{\delta}}\Delta\right)\Delta-2\mathrm{i}\alpha\frac{\partial{U}}{\partial
z}\frac{\partial}{\partial z}\right]\hat{w}^{+}(z,t)$
$\displaystyle=-\mathrm{i}\beta\frac{\partial U}{\partial
z}\hat{\eta}^{+}(z,t),$ (32)
$\displaystyle\left(\frac{2}{\Rey_{\delta}}\frac{\partial}{\partial
t}+\mathrm{i}\alpha U-\frac{1}{\Rey_{\delta}}\Delta\right)\hat{\eta}^{+}(z,t)$
$\displaystyle=0,$ (33)
$\displaystyle\hat{w}^{+}(z,T_{f})=-\frac{1}{2k^{2}}\frac{\hat{w}(z,T_{f})}{E(\mbox{\boldmath$\hat{q}$}_{o})},~{}~{}~{}\hat{\eta}^{+}(z,T_{f})$
$\displaystyle=\frac{1}{2k^{2}}\frac{\hat{\eta}(z,T_{f})}{E(\mbox{\boldmath$\hat{q}$}_{o})},$
(34) $\displaystyle\hat{w}^{+}(0,t)=\frac{\partial\hat{w}^{+}}{\partial
z}(0,t)=\hat{\eta}^{+}(0,t)$ $\displaystyle=0,$ (35)
$\displaystyle\hat{w}^{+}(z\rightarrow\infty,t)=\frac{\partial\hat{w}^{+}}{\partial
z}(z\rightarrow\infty,t)=\hat{\eta}^{+}(z\rightarrow\infty,t)$
$\displaystyle=0.$ (36)

The reader is directed to Schmid & Henningson (2001) for a thorough derivation
of equations (32) and (33). Furthermore, (31) corresponds to the expressions
to calculate the optimal forcing configuration using the adjoint fields, i.e.,

$\displaystyle\hat{f}_{u}^{opt}(z)$
$\displaystyle=-\frac{1}{2\sigma}\int_{-\infty}^{T_{f}}\left(\mathrm{i}\alpha\frac{\partial\hat{w}^{+}}{\partial
z}+\mathrm{i}\beta\hat{\eta}^{+}\right)\mathrm{e}^{-\mathrm{i}\omega_{f}t}\mathrm{d}t,$
(37) $\displaystyle\hat{f}_{v}^{opt}(z)$
$\displaystyle=\frac{1}{2\sigma}\int_{-\infty}^{T_{f}}\left(-\mathrm{i}\beta\frac{\partial\hat{w}^{+}}{\partial
z}+\mathrm{i}\alpha\hat{\eta}^{+}\right)\mathrm{e}^{-\mathrm{i}\omega_{f}t}\mathrm{d}t,$
(38) $\displaystyle\hat{f}_{w}^{opt}(z)$
$\displaystyle=-\frac{1}{2\sigma}\int_{-\infty}^{T_{f}}k^{2}\hat{w}^{+}\mathrm{e}^{-\mathrm{i}\omega_{f}t}\mathrm{d}t,$
(39)

where $\sigma$ is a Lagrange multiplier, cf. Appendix A for the details of the
derivation of (37)–(39).

The optimization problem in (29) is now transformed to a set of equations with
(15)–(18) and (21)–(25) being the state or forward equations, (32)–(36) being
the adjoint equations, and (37)–(39) being the design equations. These
equations are solved in a sequential fashion using a simple adjoint-looping
algorithm (Andersson et al., 1999). The algorithm starts with an initial guess
of $\mbox{\boldmath$\hat{f}$}^{opt}$ and iterates over the following
successive steps: (i) calculation of $\mbox{\boldmath$\hat{q}$}_{o}$ using
(15)–(18); (ii) forward-in-time integration of the state equations (21)–(25);
(ii) backward-in-time integration of the adjoint equations in (32) and (36);
(iii) updating the control terms with the available adjoint fields using
(37)–(39) and $\|\mbox{\boldmath$\hat{f}$}^{opt}\|=1$.

The forward and adjoint equations are discretized in space using a spectral
method based on Chebyshev polynomials. In this method, the equations are
mapped to the domain $\xi\in[-1,1]$ and the Gauss–Lobatto collocation
technique is utilized to obtain the discrete set of equations. This is
implemented using the differentiation matrices developed by Weideman & Reddy
(2000). Converged results are obtained for a domain size $z\in[0,20]$ and
resolution of $N_{z}=61$ Chebyshev collocation points in the vertical
direction. The initial time to start to simulations is selected to be
$t_{i}=-10\pi$. At this phase the effect of the wave is negligible. The
Crank–Nicolson scheme is employed for time integration. The time step size is
$\delta t=\pi/480$. A sensitivity analysis confirmed that the selected spatial
and temporal resolutions are sufficient.

###  3.2 Linear response of the flow

In this section, we study the linear response of the flow to the optimal
perturbations at different $\alpha$, $\beta$, $\omega_{f}$, $T_{f}$ and
$\Rey_{\delta}$. Figure 3 shows the maximum gain among all excitation
frequencies for each wavenumber pair $(\alpha,\beta)$ at a moderate Reynolds
number $\Rey_{\delta}=2000$. $\Rey_{\delta}=2000$ is the highest Reynolds
number achieved in the oscillating water tunnel of Sumer et al. (2010), where
they observed turbulent spots at the end of FPG stage. In order to study the
receptivity of SWBL among the FPG stage, the terminal time is selected to be
$T_{f}=0$. It is observed in figure 3 that SWBL is very receptive to
streamwise-constant ($\alpha=0,\beta\neq 0$) excitation in the FPG stage and
has a very weak response to two-dimensional $\alpha\neq 0,\beta=0$ and oblique
$\alpha\neq 0,\beta\neq 0$ excitations. These modes, mainly two-dimensional
ones, become only dominant in mid to late APG stage with the flow reversal.
Therefore, they do not play a role in an early subcritical bypass transition,
and will not be discussed in the remainder of the text.

Figure 3: Contours of maximum response
$G_{f}(\alpha,\beta,\omega_{f}=\omega_{f}^{m},T_{f}=0,\Rey_{\delta}=2000)$,
cf. (29), with respect to spanwise wavenumbers and terminal phase for
optimization, where $\omega_{f}^{m}$ is the excitation frequency delivering
the maximum gain.

Figure 4 further demonstrates the response of the flow to streamwise-constant
excitation on a $\beta-\omega_{f}$ plane at several terminal times
($T_{f}=-\pi/3,-\pi/6,0$, and $\pi/6$) at $\Rey_{\delta}=2000$. We see that
with increasing terminal time the frequency band to which the flow is
sensitive narrows down. Similarly, there is a shift to lower wavenumbers,
which can be linked to the growth of the boundary layer (cf. figure 2). In
general, there is a good flow response in $\beta\in[1.5,2.5]$ and
$\omega_{f}\in[0,3]$. In this range, the boundary layer amplifies the external
disturbances up to about $10^{4}$ times from the start of the wave event until
a phase at the start of APG stage ($T_{f}=\pi/6$), cf. figure 4d.

Figure 4: Contours of gain
$G_{f}(\alpha=0,\beta,\omega_{f},T_{f},\Rey_{\delta}=2000)$ with respect to
spanwise wavenumbers and excitation frequency. (a) $T_{f}=-\pi/3$; (b)
$T_{f}=-\pi/6$; (c) $T_{f}=0$; (d) $T_{f}=\pi/3$.

In order to analyze the scaling of the governing equations with Reynolds
number in the case of streamwise-constant excitation, we introduce the
transformations

$\overline{\overline{\eta}}=\frac{\hat{\eta}}{\Rey_{\delta}^{2}},~{}~{}~{}\overline{w}=\frac{\hat{w}}{\Rey_{\delta}},$
(40)

and substitude to streamwise-constant OSS equations (21) and (22), where the
terms with $\alpha$ vanish, i.e.,

$\displaystyle 2\left(\frac{\partial}{\partial
t}-\frac{1}{2}\hat{\Delta}\right)\hat{\Delta}\overline{w}$
$\displaystyle=(-\mathrm{i}\beta\frac{\partial\hat{f}_{v}}{\partial
z}-\beta^{2}\hat{f}_{w})\mathrm{e}^{\mathrm{i}\omega_{f}t},$ (41)
$\displaystyle 2\left(\frac{\partial}{\partial
t}-\frac{1}{2}\hat{\Delta}\right)\overline{\overline{\eta}}$
$\displaystyle=\frac{1}{\Rey_{\delta}}\mathrm{i}\beta\hat{f}_{u}\mathrm{e}^{\mathrm{i}\omega_{f}t}-\mathrm{i}\beta\frac{\partial
U}{\partial z}\overline{w}.$ (42)

In the scaled streamwise-constant setting, the velocity components becomes

$\overline{\overline{u}}=\frac{\mathrm{i}}{\beta}\overline{\overline{\eta}},~{}~{}~{}\overline{v}=\frac{\mathrm{i}}{\beta}\frac{\partial{\overline{w}}}{\partial
z}.$ (43)

Therefore, the evolution of cross-stream momentum by the velocity components
$v$ and $w$ is embedded into (41) and the evolution of streamwise momentum by
$u$ is linked to (42). As shown in (41) in the streamwise-constant setting the
cross-stream momentum completely decouples from the streamwise momentum. Thus,
the cross-stream perturbations are not influenced by the base flow and the
wave has no effect on them. The lack of interaction with the wave results in
the transverse components remain in their initial state, i.e.,

$\hat{v}=\hat{v}_{o}\mathrm{e}^{\mathrm{i}\omega_{f}t},~{}~{}~{}\hat{w}=\hat{w}_{o}\mathrm{e}^{\mathrm{i}\omega_{f}t}.$
(44)

Therefore, the increase in $G_{f}$ is solely due to intrinsic amplification of
$\hat{u}$ by the boundary layer.

Equation (41) suggests that introducing $\overline{w}$ rendered the cross-
stream momentum balance independent of Reynold number, while (42) indicates
that the streamwise forcing is one order lower ($O(1/Re_{\delta})$) than the
other $O(1)$ terms . Therefore, in high Reynolds numbers, direct streamwise
forcing is inefficient, and optimal external force should concentrate driving
cross-stream components, i.e.,

$\|\hat{f}_{v}^{opt}\|\approx\|\hat{f}_{w}^{opt}\|\gg\|\hat{f}_{u}^{opt}\|.$
(45)

Consequently, the streamwise forcing in (42) can be neglected for
$\Rey_{\delta}\gg 1$, and the evolution of $\overline{\overline{u}}$ becomes
independent of Reynolds number. Figure 5 validates these Reynolds-number
scalings using the numerical results for the case
$T_{f}=0,\alpha=0,\beta=1.5$, and $\omega_{f}=0$. Similar results are also
applicable to other cases. As displayed in figure 5a the streamwise component
of the optimal force is smaller than the transverse components and it vanishes
with increasing Reynolds number. Therefore, the terminal streamwise velocity
scaled with $Re_{\delta}^{2}$ and the terminal vertical velocity scaled with
$Re_{\delta}$ collapse for different Reynolds numbers, cf. figures 5b and 5c.

We now turn to input and output configurations. The optimal steady streamwise-
constant forcing configuration
$\mbox{\boldmath$f$}^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0)$ and the
resulting flow response at the terminal time are shown in the physical space
in figures 6a and 6b, respectively. In figure 6b, the contour lines present
the streamfunction defined as follows

$\tilde{v}=-\frac{\partial\tilde{\psi}}{\partial
z},~{}~{}~{}\tilde{w}=\frac{\partial\tilde{\psi}}{\partial y}.$ (46)

It is observed that the steady forcing is organized as counter-rotating cells
$(0,\hat{f}_{v}^{opt},\hat{f}_{w}^{opt})$, which induce steady rollers
$(0,\tilde{v}_{o},\tilde{w}_{o})$ with the same sense of rotation. The rollers
redistribute the streamwise momentum of the base flow, while they lift up the
low momentum fluid and pull down the high momentum fluid. As a result, streaks
that are antiphase, with the vertical velocity are produced, i.e., regions of
negative $\tilde{u}$ and positive $\tilde{w}_{o}$, and vice versa, collapse.
There is no feedback from streaks to rollers, as long as the streaks remain
stable. We will see later that the same observation also applies to nonlinear
streamwise-constant equations. Streaks are merely forced by the linear
interaction between the base flow and vertical perturbations. In (41) and
(42), $\overline{\overline{\eta}}$ and $\overline{w}$ are of the same order in
Reynolds number. So, a transverse steady forcing with an amplitude of
$O(1/Re_{\delta}^{2})$ will induce steady rollers of strength
$O(1/Re_{\delta})$. These rollers then in turn force the streaks of $O(1)$,
which will grow with increasing rates associated with the outer-velocity time-
scales. Waleffe (1995) derived a similar streak–roller system, where $O(1)$
streaks synthetically forced by the steady rollers of magnitude
$O(1/Re_{\delta})$. These scalings in Reynolds number also apply to the
streaks forced by optimal initial perturbations in the form of counter-
rotating rollers in steady boundary layers (Gustavsson, 1991; Schmid &
Henningson, 2001) and in SWBLs (Verschaeve et al., 2017).

Figure 5:  Reynolds-number dependency of the optimal linear input and output
fields. (a) Components of the optimal forcing
$\mbox{\boldmath$f$}^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0)$ at
($\cdots$): $\Rey_{\delta}=125$; ($--$): $\Rey_{\delta}=500$; ($-$):
$\Rey_{\delta}=2000$. See legends for the color coding of forcing components.
(b) Streamwise velocity $|\hat{u}(t=0)|/\Rey_{\delta}^{2}$ at the terminal
time $t=T_{f}=0$. (c) Vertical velocity $|\hat{w}_{o}|/\Rey_{\delta}$, which
is steady under steady forcing.

Figure 6:  Optimal linear input and output configurations in the physical
space. (a) Cross-stream components of the optimal steady streamwise-constant
force
$\mbox{\boldmath$f$}^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0,\Rey_{\delta}=2000)$.
Filled contours show the forcing magnitude
$|\mbox{\boldmath$f$}^{opt}|/\max\\{\mbox{\boldmath$f$}^{opt}\\}$. Streamwise
component is negligible. Arrows show
$f^{opt}_{v}\hat{\mbox{\boldmath$j$}}+f^{opt}_{w}\hat{\mbox{\boldmath$k$}}$,
where $\hat{\mbox{\boldmath$j$}}$ are $\hat{\mbox{\boldmath$k$}}$ are the
Cartesian unit vectors in spanwise and vertical directions, respectively. This
is the forcing configuration employed for the analysis in §4-6. (b) The flow
response at the terminal time $t=T_{f}=0$. Filled contours show levels of
streamwise component $\tilde{u}/\Rey_{\delta}^{2}$ and line contours show the
steady streamfunction $\tilde{\psi}_{o}/\Rey_{\delta}$ spanning 9 levels
between minimum and maximum values in the plane.

##  4 Nonlinear streaks

When the perturbations reach appreciable amplitudes, nonlinear effects should
be taken into account. We showed in § 3.2 that the cases with $\beta\neq 0$,
$\alpha=0$ and $\omega_{f}=0$ present a good balance between the optimality
and simplicity. Therefore, hereafter the discussion will focus on optimal
steady streamwise-constant perturbations, which are arranged as streaks and
rollers. In this configuration, the forcing concentrates in cross-stream
components and induces rollers that remain steady also in nonlinear regimes
due to lack of interaction with the wave. Therefore, the velocity field of the
nonlinear rollers is found by

$\displaystyle-\frac{1}{\Rey_{\delta}}\Delta\tilde{v}_{o}$
$\displaystyle=-\frac{\partial\tilde{p}_{o}}{\partial
y}+A_{o}f_{v}^{opt}-\left(\tilde{v}_{o}\frac{\partial\tilde{v}_{o}}{\partial
y}+\tilde{w}_{o}\frac{\partial\tilde{v}_{o}}{\partial z}\right),$ (47)
$\displaystyle-\frac{1}{\Rey_{\delta}}\Delta\tilde{w}_{o}$
$\displaystyle=-\frac{\partial\tilde{p}_{o}}{\partial
z}+A_{o}f_{w}^{opt}-\left(\tilde{v}_{o}\frac{\partial\tilde{w}_{o}}{\partial
y}+\tilde{w}_{o}\frac{\partial\tilde{w}_{o}}{\partial z}\right),$ (48)
$\displaystyle\frac{\partial\tilde{v}_{o}}{\partial
y}+\frac{\partial\tilde{w}_{o}}{\partial z}$ $\displaystyle=0,$ (49)

where $A_{o}$ is a small forcing magnitude with $A_{o}\ll 1$ and $\Delta$ is
the Laplacian operator. The steady rollers excite the streaks via intermodal
nonlinear interactions and also via linear interaction with the base flow,
i.e.,

$\displaystyle\frac{2}{\Rey_{\delta}}\left(\frac{\partial}{\partial
t}-\frac{1}{2}\Delta\right)\tilde{u}$
$\displaystyle=-\left(\tilde{v}_{o}\frac{\partial\tilde{u}}{\partial
y}+\tilde{w}_{o}\frac{\partial\tilde{u}}{\partial
z}+\tilde{w}_{o}\frac{\partial U}{\partial z}\right),$ (50)

where the small streamwise forcing $A_{o}f_{u}^{opt}$ is neglected. As the
evolution of the cross-stream momentum remains decoupled from the streamwise
momentum also in the nonlinear regime, there is no feedback from nonlinear
streaks to rollers as long as streaks go through streamwise-constant
deformations, which is the case for the stream-constant excitation. More
generic perturbations are to be introduced in § 5. Before proceeding with the
results for this nonlinear streak–roller system, we first transform the
nonlinear equations to a more convenient form with the aim of reducing the
number of parameters in the analysis. To this end, we introduce the variable

$A=A_{o}\Rey_{\delta}^{2},$ (51)

and define the transformations

$\accentset{\approx}{v},\accentset{\approx}{w},\accentset{\approx}{p}=\frac{\Rey_{\delta}}{A}\tilde{v},~{}\frac{\Rey_{\delta}}{A}\tilde{w},~{}\frac{\Rey_{\delta}^{2}}{A}\tilde{p}.$
(52)

Introducing the transformed variables to the roller equations

$\displaystyle-\Delta\accentset{\approx}{v}_{o}$
$\displaystyle=-\frac{\partial\accentset{\approx}{p}_{o}}{\partial
y}+f_{v}^{opt}-A\left(\accentset{\approx}{v}_{o}\frac{\partial\accentset{\approx}{v}_{o}}{\partial
y}+\accentset{\approx}{w}_{o}\frac{\partial\accentset{\approx}{v}_{o}}{\partial
z}\right),$ (53) $\displaystyle-\Delta\accentset{\approx}{w}_{o}$
$\displaystyle=-\frac{\partial\accentset{\approx}{p}_{o}}{\partial
z}+f_{w}^{opt}-A\left(\accentset{\approx}{v}_{o}\frac{\partial\accentset{\approx}{w}_{o}}{\partial
y}+\accentset{\approx}{w}_{o}\frac{\partial\accentset{\approx}{w}_{o}}{\partial
z}\right),$ (54)
$\displaystyle\frac{\partial\accentset{\approx}{v}_{o}}{\partial
y}+\frac{\partial\accentset{\approx}{w}_{o}}{\partial z}$ $\displaystyle=0,$
(55)

and to the streak equation

$\displaystyle\left(\frac{\partial}{\partial
t}-\frac{1}{2}\Delta\right)\tilde{u}$
$\displaystyle=-\frac{A}{2}\left(\accentset{\approx}{v}_{o}\frac{\partial\tilde{u}}{\partial
y}+\accentset{\approx}{w}_{o}\frac{\partial\tilde{u}}{\partial
z}+\accentset{\approx}{w}_{o}\frac{\partial U}{\partial z}\right),$ (56)

Transforming the nonlinear governing equations from (47)–(50) to (53)–(56)
reduces the parameter space of the problem from two, $\Rey_{\delta}$ and
$A_{0}$, to one, $A$, which can be considered now as the effective amplitude
of the excitation. We reiterate that this one-parameter model is only
applicable in the range of $\Rey_{\delta}\gg 1$, where optimal forcing
configuration $\mbox{\boldmath$f$}^{opt}$ does not depend on $\Rey_{\delta}$
and has a vanishing streamwise component.

The streak–roller equations are solved using the open-source CFD library
Nektar++ (Cantwell et al., 2015). To this end, a high-order spectral element
method is employed in a two-dimensional computational domain extending to
$z\in[0,L_{z}=20]$ in the vertical direction, and to $y\in[0,Ly=2\pi/\beta]$
in the spanwise direction. Periodicity is applied in the $y$ direction. The
domain discretized using a structured two-dimensional grid with $N_{y}=24$ and
$N_{z}=36$ elements in spanwise and vertical directions, respectively. The
grid is clustered towards to wall, and the expansion rate of elements in the
vertical direction is set to $1.1$. Each element is equipped with two
dimensional nodal expansion bases, which are constructed using Lagrange
polynomials that are defined on Gauss–Lobatto–Legendre points (Karniadakis &
Sherwin, 2005). A polynomial order of $N_{p}=7$ is employed. The governing
equations are projected on the polynomial basis using a continuous Galerkin
method. The resulting system of differential algebraic equations is
discretized in time using an implicit second-order scheme, cf. Vos et al.
(2011) for details. Finally, the coupled linear system of equations is
segregated using a velocity-correction scheme (Karniadakis et al., 1991).

Figure 7:  (a) Variation of the roller amplitudes (58) with respect to the
forcing amplitude A (51) in the case of linearly optimal forcing with
$\alpha=0,\beta=1.5,\omega_{f}=0$, and $T_{f}=0$. The symbols mark the values
at $A=15,50$ and $100$, for which the evolution of streak amplitudes (57) is
shown in (b). The spatial distribution corresponding to the forcing with these
amplitudes are presented in figure 8. The light horizontal lines in (b) shows
two different critical streak amplitudes that are reported in the literature
for the emergence of instabilities on steady streaks. (solid line):
$A_{s}^{c}=0.152$ by Vaughan & Zaki (2011); (dashed line): $A_{s}^{c}=0.26$ by
Andersson et al. (2001).

To keep the analysis on the evolution, stability and breakdown of nonlinear
streaks in a tractable margin, a selection has to be made for a representative
spanwise wavenumber $\beta$ and terminal time $T_{f}$. To this end, $T_{f}=0$
is a good choice to obtain strong amplification during the FPG stage.
Furthermore, we see in figure 4c,d that the wavenumber $\beta=1.5$ shows good
performance for time horizons corresponding to the strongest amplifications
($T_{f}=0,\pi/6$). Therefore, we will merely consider roller perturbations
induced by optimal forcing
$\mbox{\boldmath$f$}^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0)$ in the
current and upcoming sections. This forcing configuration was shown in figure
6a above.

It is convenient to characterize the nonlinear streaks and rollers via simple
scalar measures for their amplitudes. Following Andersson et al. (2001),
amplitude of streaks is defined as the half of the difference between maximum
and minimum perturbation velocity, i.e.,

$A_{s}(t)=\frac{1}{2}\left(\max_{y,z}\\{\tilde{u}(y,z,t)\\}-\min_{y,z}\\{\tilde{u}(y,z,t)\\}\right).$
(57)

$A_{s}$ approaches to the peak of Fourier mode ($\hat{u}$) in the linear
regime. The amplitude of steady rollers can be prescribed conveniently using
the maximum vertical velocity, i.e.,

$\tilde{w}_{\max}=\max_{y,z}\\{\tilde{w}_{o}(y,z)\\}.$ (58)

Figure 7a shows the variation of the roller magnitudes with respect to
effective forcing amplitude $A$. The amplitudes are presented in a
$\Rey_{\delta}$-independent scaling, i.e.,
$\tilde{w}_{\max}\Rey_{\delta}=A\accentset{\approx}{w}_{\max}$. We see that
the rollers are in approximately linear regime for the considered range of
forcing amplitudes $A$. Figure 7b further shows the temporal evolution of
normalized streak amplitudes $A_{s}/u_{0}(t)$ for the cases $A=15,50$ and
$100$ corresponding to roller magnitudes
$\tilde{w}_{\max}=2.8/\Rey_{\delta},9.51/\Rey_{\delta}$ and
$18.78/\Rey_{\delta}$. We see that the streaks initially grow faster than
free-stream velocity and $A_{s}/u_{0}(t)$ increases until about $t=-2$.
Subsequently, there is an equilibrium stage until about $t=-0.5$, in which
streaks and free-stream velocity grow in proportion, hence $A_{s}/u_{0}(t)$
remains approximately constant. Following this phase, the normalized streak
amplitudes increase dramatically, as steady rollers keep pumping momentum into
streaks, while free-stream velocity stagnates and decelerates. The critical
streak amplitudes calculated by Vaughan & Zaki (2011) ($A_{s}^{c}=0.152$) and
by Andersson et al. (2001) ($A_{s}^{c}=0.26$) for the initiation of
instabilities on steady streaks are also shown in figure 7b. The discrepancy
between these critical values is due to differences in the shapes of streaks
employed in these works. It is observed that the case $A=15$ remains below the
critical streak amplitudes in the FPG stage and is expected to be stable in
this stage. In contrast, cases $A=50$ and $100$ exceed the critical values
already in the FPG stage, hence can develop early instabilities. These
observations will be confirmed in § 6 using secondary stability analysis.

Figure 8:  Nonlinear streaks induced by steady streamwise-constant optimal
external excitation
$\mbox{\boldmath$f$}^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0)$ with
different amplitudes: (a) $A=0$; (b) $A=15$
($\tilde{w}_{\max}=2.8/\Rey_{\delta}$); (c) $A=50$
($\tilde{w}_{\max}=9.51/\Rey_{\delta}$); (d) $A=100$
($\tilde{w}_{\max}=18.78/\Rey_{\delta}$), cf. (51). Filled contours show
levels of total streamwise velocity scaled with the local phase value of free-
stream velocity, i.e., $U_{s}/u_{0}(t)$. Each colorbar on top shows the
contour levels in the panes below. The thick red contourlines show $u=0$.
Black line contours show the streamfunction $\tilde{\psi}_{o}$ spanning 9
levels between minimum and maximum values at the phase.

The streaky fields induced by linearly optimal steady streamwise-constant
forcing are two-dimensional and have three components, i.e.,

$[U_{s},V_{s},W_{s}](y,z,t):=[U,0,0](z,t)+[\tilde{u},0,0](y,z,t)+[0,\tilde{v}_{o},\tilde{w}_{o}](y,z),$
(59)

where the velocity components associated with the streamwise-constant
streak–roller system,
$\mbox{\boldmath$\tilde{u}$}=[\tilde{u},\tilde{v}_{o},\tilde{w}_{o}]$, added
to the standard laminar profile ($U$) of the SWBL. The spatial organization of
these fields corresponding to the cases $A=15,50$ and $100$ along with the
baseline case ($A=0$) are shown in figure 8. Filled contours in figure 8 show
the streamwise velocity fields scaled with the local free-stream velocity at
the respective phase, $U_{s}/u_{0}(t)$, and line contours show levels of
$\accentset{\approx}{\psi}_{o}$. Increasing nonlinearity with increasing $A$
leads to more uneven streak growth, where low-speed streaks are lifted up to
higher fluid layers and narrow down, e.g., compare the figures 8b,c,d. These
elevated low-momentum regions, low-speed streaks, are bounded by internal
shear layers with strong local spanwise and vertical variations.

In the case of steady streamwise-constant excitation, the gain in the
nonlinear regime reads

$G_{f}(t;\beta,A,\Rey_{\delta},\alpha=0,\omega_{f}=0)=\frac{E_{\mathcal{V}}(t)}{E_{\mathcal{V}o}},$
(60)

where

$E_{\mathcal{V}}(t)=\frac{1}{2}\int_{0}^{\infty}\int_{0}^{2\pi/\beta}[\tilde{u}^{2}+\tilde{v}_{o}^{2}+\tilde{w}_{o}^{2}]\mathrm{d}y\mathrm{d}z,~{}~{}~{}E_{\mathcal{V}o}=\frac{1}{2}\int_{0}^{\infty}\int_{0}^{2\pi/\beta}\left[\tilde{v}_{o}^{2}+\tilde{w}_{o}^{2}\right]\mathrm{d}y\mathrm{d}z,$
(61)

are the integrated kinetic energy at a time $t$ and the initial energy,
respectively. Although the streak–roller system given by
$[\tilde{u},\accentset{\approx}{v}_{o},\accentset{\approx}{w}_{o}]$ fields is
$\Rey_{\delta}$-independent, the actual flow is $\Rey_{\delta}$-dependent, as
the physical cross-stream components are
$[\tilde{v}_{o},\tilde{w}_{o}]=A/\Rey_{\delta}[\accentset{\approx}{v},\accentset{\approx}{w}]$.
Therefore, similarity with respect to $A$ only applies to streaks not to
rollers, hence the dependency of $G_{f}$ on Reynolds number. We note that
$\tilde{v}_{o}^{2}$ and $\tilde{w}_{o}^{2}$ are two orders lower in Reynolds
number than $\tilde{u}^{2}$ and therefore, can be neglected in sufficiently
high Reynolds numbers, i.e., $E_{\mathcal{V},u}\approx E_{\mathcal{V}}$, where
$E_{\mathcal{V},u}$ is the integrated streamwise kinetic energy. Consequently,
if we define a gain with similarity variables

$\tilde{G}_{f}(t;\beta,A,\alpha=0,\omega_{f}=0)=\frac{E_{\mathcal{V},u}(t)}{\tilde{E}_{\mathcal{V}o}},$
(62)

where

$\tilde{E}_{\mathcal{V}o}=\frac{1}{2}\int_{0}^{\infty}\int_{0}^{2\pi/\beta}\left[\accentset{\approx}{v}_{o}^{2}+\accentset{\approx}{w}_{o}^{2}\right]\mathrm{d}y\mathrm{d}z,$
(63)

then the quadratic dependency of $G_{f}$ on Reynolds number can be shown
explicitly

$G_{f}=\left(\frac{\Rey_{\delta}}{A}\right)^{2}\tilde{G}_{f}.$ (64)

Figure 9 shows the gains for $A=15,50$ and $100$ at $\Rey_{\delta}=2000$. The
nonlinear saturation greatly limits the amplification of streaks in higher
amplitudes, and therefore, the gains are much lower. The streaks amplify until
about $t\approx 0.9$, and then decay with the flow reversal.

Figure 9:  Nonlinear flow response measured by the gain (62) for the cases
$A=15,50$ and $100$ at $\Rey_{\delta}=2000$. The velocity fields for these
cases are presented in figures 8b-d.

##  5 Secondary instability of the nonlinear streaks

The modulated SWBL featuring the streaks presents a new laminar state, which
can be analyzed for its linear stability. In zero-pressure-gradient (ZPG)
boundary layers, once the nonlinearity saturates the streaks, their evolution
downstream is slow. Therefore, the streamwise velocity on a downstream cross
section is assumed steady and streamwise-invariant, and employed as the base
state in the secondary stability analysis (Andersson et al., 2001). We examine
the stability of nonlinear streaks in SWBL using a similar approach. The main
challenge in adapting a secondary stability analysis to the present problem is
the transient nature of SWBL – streaky base states evolve strongly under the
effect of strong dynamic and aperiodic pressure gradients. In this regard, a
suitable approach to identify temporally unstable regions beneath the wave is
the quasi-static stability analysis, in which each instantaneous state is
analyzed separately for momentary instabilities (Chen & Kirchner, 1971). Such
an assumption is only valid, if the growth rate of the instability is
significantly higher than that of the mean flow. In this regard, the quasi-
static assumption is not ideal to draw instability balloons of a transient
flow, as close to critical conditions the growth rates become too small to
satisfy the quasi-static assumption (Von Kerczek & Davis, 1974). Nevertheless,
the approach is quite useful to identify rapidly growing instabilities
relevant for transition to turbulence in the SWBL. The quasi-static assumption
for the present flow will be validated in § 6 using DNS.

Blondeaux et al. (2012) applied the quasi-static stability analysis to SWBL by
considering two-dimensional perturbations growing on one-dimensional, one-
component base profiles $U(z,t)$. Here, the analysis is extended to two-
dimensional streaky fields with three components $[U_{s},V_{s},W_{s}](y,z,t)$
as shown in (59). We consider three-dimensional tertiary perturbations of the
form

$[u^{\prime},v^{\prime},w^{\prime},p^{\prime}](x,y,z,t)=\mathrm{Re}\\{[\hat{u}^{\prime},\hat{v}^{\prime},\hat{w}^{\prime},\hat{p}^{\prime}](y,z,t)\mathrm{e}^{\mathrm{i}(\alpha
x-\int_{0}^{t}\omega(\tau)\mathrm{d}\tau)}\\},$ (65)

where $\alpha$ are real wavenumbers, and
$\omega=\omega_{\mathrm{r}}+\mathrm{i}\omega_{\mathrm{i}}$ are associated
complex frequencies. Introducing these perturbations to incompressible
Navier–Stokes equations, and linearizing the resulting equations around the
two-dimensional frozen base state $[U_{s},V_{s},W_{s}](y,z,t)$, we obtain

$\displaystyle-\frac{2}{\Rey_{\delta}}\mathrm{i}\omega\hat{u}^{\prime}+\mathrm{i}\alpha
U_{s}\hat{u}^{\prime}+V_{s}\frac{\partial\hat{u}^{\prime}}{\partial
y}+W_{s}\frac{\partial\hat{u}^{\prime}}{\partial
z}+\hat{v}^{\prime}\frac{\partial U_{s}}{\partial
y}+\hat{w}^{\prime}\frac{\partial U_{s}}{\partial
z}-\frac{1}{\Rey_{\delta}}\Delta\hat{u}^{\prime}$
$\displaystyle=-\mathrm{i}\alpha\hat{p}^{\prime},$ (66)
$\displaystyle-\frac{2}{\Rey_{\delta}}\mathrm{i}\omega\hat{v}^{\prime}+\mathrm{i}\alpha
U_{s}\hat{v}^{\prime}+V_{s}\frac{\partial\hat{v}^{\prime}}{\partial
y}+W_{s}\frac{\partial\hat{v}^{\prime}}{\partial
z}+\hat{v}^{\prime}\frac{\partial V_{s}}{\partial
y}+\hat{w}^{\prime}\frac{\partial V_{s}}{\partial
z}-\frac{1}{\Rey_{\delta}}\Delta\hat{v}^{\prime}$
$\displaystyle=-\frac{\partial\hat{p}^{\prime}}{\partial y},$ (67)
$\displaystyle-\frac{2}{\Rey_{\delta}}\mathrm{i}\omega\hat{w}^{\prime}+\mathrm{i}\alpha
U_{s}\hat{w}^{\prime}+V_{s}\frac{\partial\hat{w}^{\prime}}{\partial
y}+W_{s}\frac{\partial\hat{w}^{\prime}}{\partial
z}+\hat{v}^{\prime}\frac{\partial W_{s}}{\partial
y}+\hat{w}^{\prime}\frac{\partial W_{s}}{\partial
z}-\frac{1}{\Rey_{\delta}}\Delta\hat{w}^{\prime}$
$\displaystyle=-\frac{\partial\hat{p}^{\prime}}{\partial z},$ (68)
$\displaystyle\mathrm{i}\alpha\hat{u}^{\prime}+\frac{\partial\hat{v}^{\prime}}{\partial
y}+\frac{\partial\hat{w}^{\prime}}{\partial z}=0$ , (69)

These equations are complemented with boundary conditions
$\hat{u}^{\prime}=\hat{v}^{\prime}=\hat{w}^{\prime}=0$ on $z=0$ and
$z\rightarrow\infty$. The system of equations (66)-(69) can be written shortly
as
$L^{\prime}(U(t))\mbox{\boldmath$\hat{q}$}^{\prime}=\omega(t)\mbox{\boldmath$\hat{q}$}^{\prime}$,
where
$\mbox{\boldmath$\hat{q}$}^{\prime}(y,z,t)=[\hat{u}^{\prime},\hat{v}^{\prime},\hat{w}^{\prime},\hat{p}^{\prime}](y,z,t)$
. Using the quasi-static approximation, an eigenvalue problem is defined by
freezing the operator $L^{\prime}$ at a temporal station $t=t_{s}$ and solving
for $\omega(t_{s})$ and $\mbox{\boldmath$\hat{q}$}^{\prime}(y,z,t_{s})$, the
eigenvalue and associated eigenfunction at $t_{s}$. The solution of the
eigenproblem at each phase is obtained with Nektar++ using an Arnoldi
algorithm developed by Tuckerman & Barkley (2000) and Barkley et al. (2008).
The reader is referred to Rocco (2014) for the details. For several
representative cases, we have cross-checked the calculated leading eigenvalues
with the ones obtained by simple power iterations, where the linear equations
with random initial conditions are integrated in time until convergence to the
leading eigenvalue is achieved. Excellent agreements are found for the
imaginary parts of eigenvalues (growth rates) but noticeable differences in
real parts (frequencies) are observed. Since real parts obtained with Arnoldi
method appeared to be more sensitive to configuration parameters, we shall use
the results from power iterations in the presentation of phase velocities.

Symmetries in the two-dimensional streaky fields allow six different groups of
secondary perturbations: symmetric or antisymmetric (with respect to base-flow
streak) perturbations associated with fundamental, subharmonic and detuned
modes. Fundamental modes share the same spanwise periodicity with the base
flow, subharmonic modes have twice the periodicity of the base flow, and
detuned modes corresponds to the remaining modes, cf. Reddy et al. (1998) for
mathematical details. In ZPG boundary layers, the most growing eigenvalues are
associated with sinuous perturbations. Using an inviscid analysis, Andersson
et al. (2001) reported that fundamental and subharmonic sinuous modes have
comparable growth rates with subharmonic modes. The eigenfunctions of both
fundamental and subharmonic sinuous modes concentrate in regions on the
elevated shear layers around low-speed streaks with very similar patterns, but
fundamental modes are slightly more localized (figures 12a and 12b in
Andersson et al. (2001)). Ricco et al. (2011) employed a more comprehensive
model accounting for the effects of spatial growth and unsteadiness of
streaks, and found that fundamental sinuous modes are the most unstable modes.
In stochastic FST-driven bypass-transition scenarios, the transition is
usually initiated by breakdown of a single streak (Hack & Zaki, 2014).
Therefore, the fundamental instabilities, with their more localized nature
around the base streak, are possibly more representative for practical
situations. In this regard, we consider only fundamental-mode instabilities in
the present analysis, which allows us to use a periodic domain with a single
streak. The spatial discretization on $y-z$ plane is identical with § 4. Only
leading eigenvalues and eigenmodes are calculated.

Figure 10:  Stability of the SWBL perturbed by the linearly optimal excitation
$f^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0)$. (a) Free-stream velocity.
(b) Contours of leading imaginary eigenvalues at $\Rey_{\delta}=2000$
calculated with separate stability analysis at each $(A,t)$ using quasi-static
assumption. The presented values are the maximum values along all streamwise
wavenumbers ($\alpha$). Thick contourlines show
$\omega_{\mathrm{i}}/\Rey_{\delta}=0.0001$. (c) The maximum growth rates of
the nonlinear streaks forced with $A=0,15,50$ and $100$ (cf. figure 8). ($-$):
$\Rey_{\delta}=2000$; ($--$): $\Rey_{\delta}=4000$.

As in § 4, we only consider nonlinear streaks induced by streamwise-constant
time-invariant excitation
$f^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0,\Rey_{\delta})$. Figure 10b
shows the maximum leading eigenvalues along all streamwise wavenumbers,
$\omega_{\mathrm{i}}^{\max}(t)=\max_{\alpha}\\{\omega_{\mathrm{i}}(\alpha,t)\\}$
for varying excitation amplitudes $A$, (51). The stability boundaries are
demonstrated with thick contourlines in figure 10b. A slightly positive value
$\omega_{\mathrm{i}}^{\max}/\Rey_{\delta}=0.0001$ is employed, as the exact
neutral curve ($\omega_{\mathrm{i}}^{\max}=0$) is quite noisy in early times.
It should be stressed that the exact location of the neutral curve has little
practical bearing, as the quasi-static assumption is only physically relevant
when the instabilities evolve faster than the base flow. For weak
perturbations in the range $A\lesssim 35$, it is observed that the boundary
layer remains stable throughout the FPG stage and becomes unstable only when
the APG stage starts, i.e., the crest of the wave passes the probing location.
With further increasing forcing amplitude ($A>35$), the instabilities spread
also to the FPG stage. The growth rates for the stronger streaks in this range
rise until the flow reversal in the early APG stage and peak roughly at a
phase when maximum streak amplitudes are observed (cf. figure 9).

The maximum growth rates calculated at $\Rey_{\delta}=2000$ and
$\Rey_{\delta}=4000$ are shown separately in figure 10c for the cases
$A=0,15,50$ and $100$, for which the streak amplitudes are plotted in figure
7b and the base states are presented in figure 8. The results for the
unperturbed case $A=0$ corresponds to the growth rates of the primary two-
dimensional instabilities. These orderly instabilities take place in the APG
stage following the emergence of inflectional profiles. The details of the
primary instabilities are well documented elsewhere, e.g., Blondeaux et al.
(2012) and Sadek et al. (2015). The case $A=15$ represents a case in the
regime with weak streak amplitudes (figure 7b). Interestingly, the secondary
instability in this case has lower growth rates than the baseline primary
instability. In the peaking phase ($t\approx 1$), the growth rate in $A=15$ is
about the half of the one in $A=0$. These reduced growth rates suggest a
damping mechanism introduced to the flow by weak-amplitude streaks. This will
be elaborated later. The last two cases, $A=50$ and $100$, showcase the
results for strong streaks that can develop instabilities already in the FPG
stage. As the FPG stage is linearly stable in the unperturbed SWBL, these
early instabilities in the FPG stage imply a possibility for a subcritical
transition. The seeding phase of the instabilities in $A=50$ and $A=100$
roughly corresponds to the phases when the streak amplitudes exceeds the
critical threshold given by Vaughan & Zaki (2011), cf. figure 7b. We further
see in figure 10c that the scaled growth rates,
$\omega_{\mathrm{i}}^{\max}/\Rey_{\delta}$, for the instabilities in $A=50$
and $100$ are independent of Reynolds number for the selected range (cf. blue
and green lines). There is a very weak dependence on Reynolds number for the
other cases, cf. the discrepancy at the start of APG stage for magenta and red
lines. The reader is referred to Sadek (2015) for the analysis of Reynolds-
number dependency of the primary two-dimensional instabilities ($A=0$).

Figure 11:  The overall growth measured by the modal perturbation kinetic
energy density $E_{\alpha}$ at $t=\pi$, cf. (70). The initial energy density
is set to $E_{0}=1$. The most unstable streamwise wavenumber ($\alpha$) is
employed at each $A$. The symbols refer to the cases $A=0,15,50$ and $100$,
which are further elaborated in figures 12 and 13.

Figure 12:  The variation of leading imaginary eigenvalues with respect to
streamwise wavenumber for the streaks induced by linearly optimal excitation
$f^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0,\Rey_{\delta}=2000)$. (a)
Baseline case with $A=0$ (51) corresponding to the roller magnitude
$\tilde{w}_{\max}=0$ (58); (b) $A=15,\tilde{w}_{\max}=0.0014$; (c)
$A=50,\tilde{w}_{\max}=0.0048$; (d) $A=100,\tilde{w}_{\max}=0.0094$. Color
scales of the contours are different in each pane and are shown in separate
colorbars next to the panes.

It was observed in figure 10 that the streaks have a dual role in the
transition to turbulence beneath solitary waves: they can be stabilizing or
destabilizing depending on their amplitude. This can be elaborated by
considering the overall growth of the perturbation energy in time. Using the
ansatz (65), the energy at a mode is expressed by

$E_{\alpha}(t)=E_{0}\mathrm{e}^{2\int_{0}^{t}\omega_{\mathrm{i}}(\tau)\mathrm{d}\tau},$
(70)

where $E_{0}$ is the initial energy density at the mode. Figure 11
demonstrates the variation of $E_{\alpha}$ at $t=\pi$ with respect to $A$. In
this figure, the most unstable $\alpha$ for each respective $A$ is evaluated.
Three different instability regimes are observed: (i) primary instability
($A=0$); (ii) inner-instability regime (say $0<A<20$, or
$0<\tilde{w}_{\max}<3.76/\Rey_{\delta}$ using figure 7a), and (iii) outer-
instability regime ($A>20$, or $\tilde{w}_{\max}>3.76/\Rey_{\delta}$). We
employed here the naming convention proposed by Vaughan & Zaki (2011), where
“inner” and “outer” refer to the location of the critical layer. In the inner-
instability regime, streaks are weak but still effective in mixing the
momentum of the base profiles and introducing a damping effect. Consequently,
there is a reduction in the growth rates of instabilities. A similar
stabilizing effect is observed in flat-plate boundary layers when moderate-
amplitude streaks are superposed on TS waves (Cossu & Brandt, 2004; Liu et
al., 2008b). The temporally unstable phases in the inner-instability regime
roughly overlaps with the baseline instabilities in the undisturbed regime
(cf. figure 10b), which suggests a modified instability of similar nature,
where the primary driving mechanism is vertical shear ($\partial
U_{s}/\partial z$). We will show later the inner instabilities develop in
regions nearer to the wall, where the vertical shear is strong, hence the name
“inner”. In the third regime, the streaks are strong enough to develop
secondary shear layer instabilities in the elevated outer zones. The overall
growth due to these instabilities rise dramatically between $20<A<60$.
Afterwards, there is a saturation range until $A\approx 90$, in which
increased forcing does not lead to substantial growth. However, with further
increasing $A$ there is another receptive regime for $A>90$.

Figure 13:  Filled red contours show four different levels of the modulus of
the streamwise component of the leading eigenmodes, $|\hat{u}^{\prime}|$,
calculated at four different times $t=-\pi/6,0,\pi/6,\pi/3$. The base states
are nonlinear streaks induced by the linearly optimal excitation
$f^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0,\Rey_{\delta}=2000)$ for
excitation magnitudes (a-d): $A=0,\alpha=0.45$; (e-h): $A=15,\alpha=0.35$;
(i-l): $A=50,\alpha=0.75$; (m-p): $A=100,\alpha=0.75$. Contour lines show the
streamwise component of the base flow scaled with the free-stream velocity at
the respective phase,
$U_{s}(y,z,\theta_{s}=-40^{\circ})/u_{0}(t)=\\{0.1:0.1:0.9\\}$. The dashed
lines represent the critical velocity level $U_{s}^{c}$, at which
$U_{s}=c_{r}$.

Figure 12 demonstrates the variation of growth rates with respect to
streamwise wavenumbers $\alpha$ for the cases $A=0,15,50$ and $100$ calculated
at $\Rey_{\delta}=2000$. It is observed that the primary instabilities have
relatively longer wavelength peaking at wavenumbers around $\alpha\approx
0.45$ (figure 12a), where the outer instabilities concentrate in shorter
wavelengths, e.g., $\alpha\approx\beta/2\approx 0.75$ for $A=50$ (figure 12c)
and $\alpha\approx\beta\approx 1.5$ for $A=100$ (figure 12d). We see that the
case $A=100$ is sensitive to a wider range of instabilities, e.g., $A=50$ is
mostly stable for the short wave perturbations with $\alpha>1.5$ whereas
$A=100$ is still very unstable at this range. It is this additional support
for highly unstable short-wavelength modes, which gives rise to a second
receptive regime for $A>90$ in figure 11. The case $A=15$, belonging to inner-
instability regime, shows a mixed behavior. Right after the flow reversal,
there is a peak at $t\approx 1$ and $\alpha\approx 0.35$ but then the growth
of the long waves stagnate (figure 12b). Meanwhile, we observe another growth
region concentrated at shorter wavelengths close to $\alpha\approx 0.75$. This
is the typical wavenumber for the short wave outer instabilities, which
suggests that outer instabilities are influential at mid to late APG stage in
the case $A=15$. However, as the overall growth ($E_{\alpha}$) at
$\alpha=0.35$ is higher, the inner-instabilities are the dominant secondary
mechanism at $A=15$.

We now turn to the the nature of instabilities, e.g., the symmetry patterns,
phase velocities, amplification mechanisms. To this end, we select the cases
with maximum growth rates in each instability regime, i.e.,
$(A=0,\alpha=0.45)$; $(A=15,\alpha=0.35)$; $(A=50,\alpha=0.75)$;
$(A=100,\alpha=0.75)$. Figure 13 demonstrates the spatial distribution of
eigenmodes using the modulus of streamwise components. It is observed that the
unstable modes in primary-instability ($A=0$) and inner-instability ($A=15$)
regimes extend to the whole spanwise extent of the periodic domain, cf.
figures 13c,d,g,h. In contrast, the eigenmodes are located around the elevated
low-speed streaks and are of more localized nature in the outer-instabilities
regime, cf. figures 13i-p. The instabilities in streaky flows are generated by
inviscid mechanisms due to inflection points in the shear layers. In these
instabilities, critical layers, where
$U_{s}^{c}:U_{s}=c_{r}=\omega_{\mathrm{r}}/\alpha$, form. The eigenmodes
concentrate in the critical layers, where they convect with local mean
velocity $U_{s}^{c}$, cf. dashed lines in figure 13. Streaky boundary layers
consist of spanwise and vertical shear layers associated with $\partial
U_{s}/\partial y$ and $\partial U_{s}/\partial z$, respectively. These are
shown in figure 14 for two representative cases $(A=15,\alpha=0.35)$ and
$(A=50,\alpha=0.75)$ at time $t=\pi/6$. We see that the critical layer in the
inner-instability regime develops on the vertical shear layer close to the
wall (figure 14b), while the critical layer in the outer-instability regime is
located in the elevated spanwise shear layers around the low-speed streak
(figure 14c).

The spanwise and vertical shear layers in streaky boundary layers have a
dampening effect, when the principle instability does not develop on them.
This can be understood by breaking down the generation of the modal
perturbation kinetic energy into individual components. Following Cossu &
Brandt (2004), we express the globally integrated budget as follows

$\frac{\partial e_{\mathcal{V}}}{\partial
t}=p_{\mathcal{V},y}+p_{\mathcal{V},z}-d_{\mathcal{V}},$ (71)

where $e_{\mathcal{V}}$ is the total perturbation kinetic energy,
$p_{\mathcal{V},y}$ is the total production rate due to spanwise shear,
$p_{\mathcal{V},w}$ is the total production rate due to vertical shear and
$d_{\mathcal{V}}$ is the total dissipation rate. Here we neglect the
contributions related to base spanwise ($V_{s}$) and vertical ($W_{s}$)
velocities as $\tilde{v}_{o}$ and $\tilde{w}_{o}$ are an order of magnitude
smaller in $\Rey_{\delta}$. If we consider the perturbations associated with a
single instability mode, we can express individual terms by
$[e_{\mathcal{V}},p_{\mathcal{V},y},p_{\mathcal{V},z},d_{\mathcal{V}}]=[\hat{e}_{\mathcal{V}},\hat{p}_{\mathcal{V},y},\hat{p}_{\mathcal{V},z},\hat{d}_{\mathcal{V}}]\mathrm{e}^{2\omega_{\mathrm{i}}t}$,
where

$[\hat{e}_{\mathcal{V}},\hat{p}_{\mathcal{V},y},\hat{p}_{\mathcal{V},z},\hat{d}_{\mathcal{V}}]=\frac{1}{\lambda_{y}}\int_{0}^{\lambda_{y}}\int_{0}^{\infty}[\hat{E},\mathcal{\hat{P}}_{v},\mathcal{\hat{P}}_{w},\hat{\mathcal{D}}]\mathrm{d}y\mathrm{d}z,$
(72)

$\lambda_{y}=2\pi/\beta$ and

$\displaystyle\begin{split}\hat{E}=\frac{2}{\Rey_{\delta}}(\hat{u}^{\prime*}\hat{u}^{\prime}+\hat{v}^{\prime*}\hat{v}^{\prime}+\hat{w}^{\prime*}\hat{w}^{\prime}),&~{}~{}~{}~{}~{}\mathcal{\hat{D}}=\frac{2}{\Rey_{\delta}}(\hat{\varepsilon}^{\prime*}\hat{\varepsilon}^{\prime}+\hat{\zeta}^{\prime*}\hat{\zeta}^{\prime}+\hat{\eta}^{\prime*}\hat{\eta}^{\prime}),\\\
\hat{P}_{y}=-(\hat{u}^{\prime*}\hat{v}^{\prime}+\hat{v}^{\prime*}\hat{u}^{\prime})\frac{\partial
U_{s}}{\partial
y},&~{}~{}~{}~{}~{}~{}~{}\hat{P}_{z}=-(\hat{u}^{\prime*}\hat{w}^{\prime}+\hat{w}^{\prime*}\hat{u}^{\prime})\frac{\partial
U_{s}}{\partial z}.\end{split}$ (73)

In the dissipation term $\hat{\mathcal{D}}$, the components of the
perturbation vorticity $[\varepsilon^{\prime},\zeta^{\prime},\eta^{\prime}]$
are employed. Now the growth rate of the instability can be expressed as

$\omega_{\mathrm{i}}=\frac{1}{2\hat{e}_{\mathcal{V}}}(\hat{p}_{\mathcal{V},y}+\hat{p}_{\mathcal{V},z}-\hat{d}_{\mathcal{V}}).$
(74)

Figure 14:  Derivative fields of the base streamwise velocity $U_{s}$ at
$t=\pi/6$. (a) $\partial U_{s}/\partial y$ in $A=15,\alpha=0.35$, (b)
$\partial U_{s}/\partial z$ in $A=15,\alpha=0.35$, (c) $\partial
U_{s}/\partial y$ in $A=50,\alpha=0.75$, (d) $\partial U_{s}/\partial z$ in
$A=50,\alpha=0.75$. The dashed lines show the critical layers, where
$U_{s}=c_{r}$. The streaks are induced by steady streamwise-constant optimal
external excitation
$\mbox{\boldmath$f$}^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0,\Rey_{\delta}=2000)$.

Figure 15:  Production and dissipation rates of perturbation kinetic energy
(cf. 73) for cases $A=15,\alpha=0.35$ (a-c) and $A=50,\alpha=0.75$ (d-f) at
time $t=\pi/6$, cf. figures 13g,k. The energy of the perturbations is
normalized to unity, i.e., $e_{\mathcal{V}}=1$, in the calculation of fields.
The streaky base states are induced by steady streamwise-constant optimal
external excitation
$\mbox{\boldmath$f$}^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0,\Rey_{\delta}=2000)$.
(a,d): Production rate due to spanwise shear ($\hat{\mathcal{P}}_{y}$); (b,e):
Production rate due to vertical shear ($\hat{\mathcal{P}}_{z}$); (c,f):
Dissipation rate ($\hat{\mathcal{D}}$). The integrated contributions of the
fields, $\hat{p}_{\mathcal{V},y},\hat{p}_{\mathcal{V},y}$ and
$\hat{d}_{\mathcal{V}}$, are also presented in the respective panels.

Figure 15 shows the production and dissipation fields for the cases
$A=15,\alpha=0.35$ and $A=50,\alpha=0.75$ at time $t=\pi/6$. The total
integrated values,
$\hat{p}_{\mathcal{V},y},\hat{p}_{\mathcal{V},y},\hat{d}_{\mathcal{V}}$, are
also presented in the respective panels of the fields. The energy of the
perturbations is normalized to unity, i.e., $e_{\mathcal{V}}=1$. In the case
of inner instability ($A=15$), the production due to vertical shear feeds the
growth (figure 15b), while the production due to spanwise shear has negative
contributions to the total budget, i.e., has a stabilizing effect (figure
15c). Vice versa is true for the outer instability – the spanwise shear drives
the instability, while vertical shear trying to counteract it, cf. figures
15d,e. The degree of dualism between the two shear production mechanisms and
the dissipation rate determines together the growth rate of the instability,
cf. (74). In case $A=15$, the shear-damping effect is stronger with
$\left|\hat{p}_{\mathcal{V},y}\right|$ being about $20\%$ of
$p_{\mathcal{V},y}$. Furthermore, the dissipation rate is also higher in this
case, as the perturbations are closer to the wall where viscous effects are
more pronounced.

The counteracting role of vertical and spanwise shear layers in boundary
layers has been well documented in steady flows. Reddy et al. (1998) discussed
the stabilizing effect of the vertical shear on the outer instabilities
developing on high-amplitude streaks. The same shear damping applies in outer-
instability regime in SWBLs (figures 15d,e). For low-amplitude streaks, Cossu
& Brandt (2004) has reported inner TS-like instabilities with reduced growth
rates. They suggested the negative contributions from
$\hat{p}_{\mathcal{V},y}$ as the primary mechanism behind the stabilizing
effect of low-amplitude streaks. This term vanishes in the unperturbed
boundary layer ($\partial U/\partial y=0$), while the vertical production
rates remain at similar magnitude, hence the higher growth rate of the
undisturbed TS instability. The same stabilization mechanism is also effective
in the inner-instability regime of SWBLs. However, we note that another
mechanism contributing to the reduction of growth rates is the increase in
dissipation due to three dimensionality. Orderly instability modes are two-
dimensional and have one component vorticity ($\zeta^{\prime}$), which yields
lower dissipation rates, cf. $\hat{\mathcal{D}}$ in (73).

Figure 16:  Phase velocities $c_{r}=\omega_{\mathrm{i}}/\alpha$ for the cases
presented in figure 13. (a) Absolute values. (b) Normalized values with the
local free-stream velocity at the respective phases. Dashed line shows
$c_{r}/u_{0}(t)=0.82$, which is the phase velocity calculated by Andersson et
al. (2001) for a outer streak instability in a ZPG boundary layer. Dotted line
shows $c_{r}/u_{0}(t)=0.31$, which is the phase velocity calculated by Cossu &
Brandt (2004) for an inner (modified TS-like) instability.

Figure 17:  Symmetry patterns of instabilities are shown using the isosurfaces
of the streamwise component of eigenmodes. (a) two-dimensional mode at
streamwise wavenumber $\alpha=0.45$ growing on baseline flow with $A=0$ at
time $t=\pi/3$ (cf. figure 13d); (b) varicose mode at $\alpha=0.35$ growing on
low-amplitude streaks with $A=15$ at $t=\pi/3$ (cf. figure 13h); (c) sinuous
mode at $\alpha=0.75$ growing on streaks with $A=50$ at $t=\pi/6$ (cf. figure
13k). (light): negative isosurface; (dark): positive isosurface.

The characteristics of the instabilities can be further elaborated by studying
the phase velocities and symmetry patterns. Figure 16 plots the phase
velocities for the cases presented in figure 13, where figure 16a is scaled
with $u_{0,m}^{*}$ and figure 16b is in local scaling with the free-stream
velocity at the respective phase ($c_{r}/u_{0}(t)$). Additionally, the phase
velocity of the outer streak instability in Andersson et al. (2001)
($c_{r}=0.82$) and the inner TS-like instability in Cossu & Brandt (2004)
($c_{r}=0.31$) are also shown in figure 16b for reference. We observe in
figure 16b that the phase velocity of outer instabilities in FPG stage has
very close values to their counterparts in ZPG boundary layers, cf. light-
colored lines in $t<0$ in figure 16b. The deceleration in the APG stage has a
dramatic effect on the phase velocity of these modes, and they decay rapidly
in the APG stage. The instabilities generated by the vertical shear ($A=0,15$)
follows initially the phase velocity calculated by Cossu & Brandt (2004) but
then decays also rapidly with flow reversal in the vicinity of the wall.

Figure 17 shows the symmetry pattern of the primary, inner and outer
instabilities using the streamwise velocity of instances shown in figures
13d,h,k. The primary instability is as expected two dimensional (cf. figure
17a). On the other hand, the inner instability is a varicose mode strongly
tilted in the streamwise direction (cf. figure 17b). A varicose symmetry is
also reported in Cossu & Brandt (2004)) for the inner TS-like instability.
Finally, we see in figure 17c that the outer instability is in the form of
sinuous mode. Sinuous modes are the most unstable mode in ZPG boundary layers
(Andersson et al., 2001; Ricco et al., 2011) and also commonly observed as
breaking mode of streaks in experiments (e.g. Mans et al. (2007)).

##  6 Direct numerical simulations

In § 4, we analyzed the dynamics of streamwise-constant streaks in a two-
dimensional domain ($\mbox{\boldmath$x$}=[x=0,y,z]$). In this section, we
investigate the response of the streaks to small-amplitude background noise to
validate the quasi-static assumption employed in § 5 and to investigate
breakdown stage in transition. To this end, three-dimensional perturbations
are introduced and direct numerical simulations are conducted. The optimal
forcing configuration remain identical to the one employed in § 4 and § 5,
i.e., $f^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0,\Rey_{\delta})$.

Case | $A$ | $A_{r}$ | $\Rey_{\delta}$ | $T_{r}$ | $\alpha_{0}$ | $\beta_{0}$ | $L_{x}\times L_{y}\times L_{z}$ | $N_{x}\times N_{y}\times N_{z}$
---|---|---|---|---|---|---|---|---
A0c1 | 0 | $10^{-1}$ | 4000 | 0 | $0.225$ | - | $2\pi/\alpha_{0}\times 4\pi/\alpha_{0}\times 20$ | $400\times 320\times 480$
A0c2 | 0 | $10^{-6}$ | 4000 | 0 | $0.225$ | - | $2\pi/\alpha_{0}\times 4\pi/\alpha_{0}\times 20$ | $400\times 320\times 480$
A0c3 | 0 | $10^{-12}$ | 4000 | 0 | $0.225$ | - | $2\pi/\alpha_{0}\times 4\pi/\alpha_{0}\times 20$ | $400\times 320\times 480$
A15c1 | 15 | $10^{-1}$ | 4000 | $0$ | $0.175$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $400\times 320\times 480$
A15c2 | 15 | $10^{-6}$ | 4000 | $0$ | $0.175$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $400\times 320\times 480$
A15c3 | 15 | $10^{-12}$ | 4000 | $0$ | $0.175$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $400\times 320\times 480$
A50c1 | 50 | $10^{-6}$ | 2000 | $-\pi$ | $0.375$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $160\times 320\times 480$
A50c2 | 50 | $10^{-12}$ | 2000 | $-\pi$ | $0.375$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $160\times 320\times 480$
A100c1 | 100 | $10^{-6}$ | 2000 | $-\pi$ | $0.375$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $160\times 320\times 480$
A100c2 | 100 | $10^{-9}$ | 2000 | $-\pi$ | $0.375$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $160\times 320\times 480$
A100c3 | 100 | $10^{-11}$ | 3000 | $-\pi$ | $0.375$ | 1.5 | $2\pi/\alpha_{0}\times 2\pi/\beta_{0}\times 20$ | $160\times 320\times 480$
Table 1: Computational details of simulations. $\beta_{0}$ is the spanwise
wavenumber of the excitation. $A_{r}$ is the amplitude of random tertiary
perturbations seeded at time $T_{r}$. Each element has in total $(P+1)^{2}$
degrees of freedom (DOF), where $P=7$ is the order of the Lagrange
polynomials. The total number of DOFs in $y$ and $z$ directions is calculated
as, e.g. for $y$ direction, $N_{y}=N_{ey}\times(P+1)$, where $N_{ey}$ is the
number of elements in the $y$ direction. In the streamwise direction,
$N_{x}/2$ Fourier modes are employed yielding $N_{x}$ grid points.

The computational details of the cases are presented in Table 1. We consider
four representative forcing amplitudes $A=0,15,50$ and $100$ as in previous
sections and perturb the resulting streaky fields with small-amplitude random
noise of amplitude of $A_{r}$. These random perturbations are seeded at the
end of FPG stage $T_{r}=0$ for the cases with $A=0$ and $A=15$, as these cases
are stable in the FPG stage. For the rest, the tertiary random perturbations
are seeded at $T_{r}=-\pi$. The numerical method employed in § 4 is extended
to three dimensions using a mixed spatial discretization in Nektar++. In this
method, a bi-dimensional spectral-element method, previously introduced in §
4, is employed in streamwise-wall normal ($y-z$) plane, and global Fourier
expansions are considered in the streamwise ($x$) direction (Karniadakis,
1990; Bolis, 2013). To avoid instabilities due to aliasing errors, the method
developed by Kirby & Sherwin (2006) is employed for polynomial expansions, and
the $2/3$ rule is employed for the Fourier expansions (Boyd, 2001). The
computational domain is a rectangular box with dimensions
$[0,L_{x}]\times[0,L_{y}]\times[0,L_{z}]$. Periodic boundary conditions are
employed in the streamwise and spanwise directions. The domain contains a
single streak in the forced cases. The streamwise extent of the domains is
selected to allow growth in the most unstable streamwise wavenumbers. No-slip
boundary condition is applied at the bottom wall, and free-slip boundary
condition is applied at the top wall. Verschaeve & Pedersen (2014) remarked
that a very fine grid resolution is necessary to capture the natural
development of two-dimensional instabilities. Therefore, a finer structured
grid than the one in § 4 is employed to resolve instabilities and turbulence.
The grid densities are everywhere considerably higher than the previous DNS
works on SWBL (Vittori & Blondeaux, 2008; Ozdemir et al., 2013).

Figure 18:  Growth and nonlinear saturation of secondary instabilities. Lines
show the modal energy extracted from DNS, whereas symbols show the ones
calculated using the leading eigenvalues of secondary stability analysis (§
5). (a) Cases A0c1, A0c2, A0c3; (b) Cases A15c1, A15c2, A15c3; (c) Cases
A50c1, A50c2; (d) Cases A100c1, A100c2, A100c3, cf. Table 1 for case
definitions. The horizontal line shows $E_{\alpha}=10^{-2}$.

The energy density in each streamwise mode $\alpha$ is calculated by Fourier
transforming the velocity fields in the streamwise direction and integrating
the respective energy in the Fourier mode
$\mbox{\boldmath$\hat{u}$}(\alpha,y,z,t)$ over the domain and then normalizing
it, i.e.,

$E_{\alpha}(t)=\frac{1}{2L_{y}L_{z}}\int_{0}^{L_{z}}\int_{0}^{L_{y}}\mbox{\boldmath$\hat{u}$}^{*}(\alpha,y,z,t)\cdot\mbox{\boldmath$\hat{u}$}(\alpha,y,z,t)\mathrm{d}y\mathrm{d}z.$
(75)

Since the introduced random perturbations are of small magnitude, we expect
linear mechanisms will drive the initial growth of secondary instabilities.
Therefore, the modal kinetic energies extracted from the direct numerical
simulations should match the ones calculated with the secondary stability
analysis with (61), if the quasi-static assumption is valid. This is tested in
figure 18 for all cases, where the initial modal energy level ($E_{0}$) of the
linear growth results is adjusted to match the DNS values. An interesting
first observation is that the energy growth in all considered cases saturates
at about $E^{c}:=E_{\alpha}=10^{-2}$ regardless of Reynolds number, saturation
phase and the type of instability. This energy level can be considered as the
critical threshold for the onset of breakdown. The cases, which cannot reach
this level during the wave event, can be assumed still laminar. We observe a
good match between DNS and linear stability theory (LST) in the cases with
outer instabilities (figure 18c,d). In these cases, the long wave instability
at $\alpha=0.375$ develops first thanks to its higher growth rate in the FPG
stage. Depending on the initial noise amplitude, these long-wave modes can
reach the critical level and become the mode of breakdown as in the Case
A50c1, or are overtaken by the shorter wave instability at $\alpha=0.75$ as in
the Case A50c2, cf. figure 18c.

There is also good agreement between DNS and LST in the cases containing
inner-instabilities (cf. figure 18b). However, the DNS data stagnates in Case
A15c3 in the interval $t>1.5$ and does not follow the growth dictated by LST
anymore (only until $t=1.9$ shown in the figure). This deviation suggests that
the instabilities introduce non-negligible deformations to the slow base flow
in the late APG stage. Thus, the quasi-static assumption appears to be
inapplicable to later phases of the wave event. The biggest discrepancy
between DNS and LST is observed in two-dimensional baseline instabilities, cf.
figure 18a. In these cases, the instabilities in DNS develop with some delay
compared to the theoretical predictions. The stabilizing effect of weak
streaks can be clearly seen in figures 18a,b. The onset of transition is
substantially delayed in cases with $A=15$ compared to those with $A=0$. In
fact, in the cases with lowest initial noise, Case A15c3 remains laminar,
whereas Case A0c3 breaks into turbulence at about $t\approx 1.3$.

If we assume that LST results are always applicable and all instabilities are
of inviscid nature with constant $\omega_{\mathrm{i}}/\Rey_{\delta}$, then we
can utilize the empirical threshold $E^{c}$ and the growth rates
$\omega_{\mathrm{i}}/\Rey_{\delta}$ calculated at a specific Reynolds number
(e.g. $\Rey_{\delta}=2000$) to extrapolate our results to a wider range of
Reynolds numbers and perturbation levels. Using this extrapolation, we can
generate state diagrams showing whether the flow is laminar or turbulent at an
instant $t$. To this end, the state of the flow is a function of four
parameters, $t$, $\tilde{w}_{\max}$, $\Rey_{\delta}$, and the initial
perturbation energy in the instability mode, $E_{0}$. Figures 19 show the flow
states with respect to $\tilde{w}_{\max}$ and $\Rey_{\delta}$ at $t=2/9\pi$
and $t=\pi/2$ for initial perturbation levels of $E_{0}=10^{-20}$ and
$E_{0}=10^{-32}$. The cases sharing the same initial perturbation levels are
also demonstrated with symbols in the respective diagrams. The boundary
between inner and outer instabilities, i.e., $A=20$ corresponding to
$\tilde{w}_{\max}=3.76/\Rey_{\delta}$ (figure 7a), is also plotted in figures.
As shown in figures 19a,b the primary and inner streak instabilities are not
effective yet at $t=2/9\pi$. At this earlier phase, the transition occurs only
due to outer instabilities, which develop when rollers exceed a certain
threshold depending on Reynolds number. This is the manifestation of bypass
transition. The primary instability modes are bypassed by an early subcritical
transition mechanism that is dependent on the magnitude of environment
perturbations, $\tilde{w}_{\max}$ in our model. The flow states are also
somewhat sensitive to the amplitude of initial tertiary perturbations
($E_{0}$) especially in lower Reynolds number, e.g., compare the range
$1000<\Rey_{\delta}<1500$ in figures 19a,b.

Figure 19:  State of the SWBL with respect to Reynolds number
($\Rey_{\delta}$) and amplitude of steady-roller perturbations
($\tilde{w}_{\max}$) are shown at two representative times in the APG stage
for two different initial tertiary perturbations ($E_{0}$). The rollers are
induced by steady streamwise-constant forcing
$f^{opt}(\alpha=0,\beta=1.5,\omega_{f}=0,T_{f}=0)$. (a) $t=2/9\pi$,
$E_{0}=10^{-32}$; (b) $t=2/9\pi$, $E_{0}=10^{-20}$; (c) $t=\pi/2$,
$E_{0}=10^{-32}$; (d) $t=\pi/2$, $E_{0}=10^{-20}$. The red dashed lines
demonstrate the boundary ($A=20$) between inner and outer instabilities, which
corresponds to $\tilde{w}_{\max}=3.76/\Rey_{\delta}$ (figure 7a).

When the wave propagates further, the primary and inner instabilities become
active, cf. figures 19c,d. We see that the laminar region protrudes into the
turbulent region in the range $\tilde{w}_{\max}\approx 0.001-0.002$, i.e., the
flow remains laminar until relatively high Reynolds numbers in this range.
This is the manifestation of the stabilization introduced by weak streaks. For
instance, Case A1c3, which has a roller magnitude of
$\tilde{w}_{\max}=2.8/\Rey_{\delta}$, remains in the protruded laminar region
for $E_{0}=10^{-32}$, cf. figure 19c. We further observe in figures 19c,d that
the primary and inner instabilities are more sensitive to $E_{0}$ compared
with outer instabilities. These results show that transition to turbulence in
the SWBL depends on the amplitude of environment perturbations even in the
case of orderly transition with two-dimensional instability modes. Flow-state
classifications that are based merely on Reynolds number, e.g., the state
charts in Sumer et al. (2010) and Ozdemir et al. (2013), have to be extended
to include some measure for environment perturbations.

Figure 20:  Streak breakdown and onset of turbulence in Case A15c1. White
surfaces show the low-speed streak using an isosurface of streamwise
fluctuation velocity $\tilde{u}^{\prime}=u-\langle u\rangle$, where $\langle
u\rangle(z,t)$ is the average value on a plane at $z$. Colored isosurfaces
show the vortical regions using Q criterion. (a) $t=23/90\pi$,
$\tilde{u}^{\prime}=-0.18$; $Q=0.003$ (b) $t=24/90\pi$,
$\tilde{u}^{\prime}=-0.18$; $Q=0.15$ ; (c) $t=25/90\pi$,
$\tilde{u}^{\prime}=-0.18$; $Q=0.55$.  Figure 21:  Streak breakdown and onset
of turbulence in Case A50c2, cf. figure 20 for the definition of surfaces. (a)
$t=50/180\pi$, $\tilde{u}^{\prime}=-0.26$; $Q=0.006$ (b) $t=52/180\pi$,
$\tilde{u}^{\prime}=-0.16$; $Q=0.027$ ; (c) $t=55/180\pi$,
$\tilde{u}^{\prime}=-0.13$; $Q=0.44$.

Figure 20 shows the breakdown of inner instability in Case A15c1. At
$t=24/90\pi$ in figure 20a, we see at the center a low-speed streak making
undulations in the downstream direction with wavenumber $\alpha=0.35$. Since
the inner instability is of varicose nature, the undulations are symmetric
with respect to the streak. Vortical structures around the low-speed streak
are also shown in the figure using a positive isosurface of Q-criterion (Hunt
et al., 1988). Among several vortical features, $\Lambda$-like vortices can
bee seen to accompany the undulating streak, cf. e.g., the region $10<x<20$ in
figure 20a. These features are reminiscent of $\Lambda$ vortices developing on
streak-modulated instability waves in ZPG boundary layers (Liu et al., 2008a).
Later at $t=25/90\pi$, the breakdown to small scales is initiated in the near-
wall layers, while the low-speed streak remains still stable and coherent, cf.
small-scale vortices in figure 20b. Subsequently, chaotic small-scale motions
quickly spread everywhere, the streak is disintegrated and the onset of
turbulence is completed at $t=26/90\pi$, cf. figure 20c.

The transition to turbulence in Case C50c2 is demonstrated in figure 21.
Initially at $t=50/180\pi$, we see a low-speed streak at the center of domain
occupying the whole streamwise extent, cf. figure 21a. This streak is unstable
and exhibits sinuous undulations with a streamwise wavelength corresponding to
the dominant outer instability mode at $\alpha=0.75$. Subsequently, at
$t=52/180\pi$, the waviness of streaks is increased and some more tertiary
vortical features have emerged, cf. figure 21b. Both vortex and velocity
structures appear to be large-scale organized features, thus the flow is still
at a laminar transitional state at this phase. Finally, at $t=55/180\pi$,
turbulence sets in and chaotic motions are to be seen everywhere in the
domain, cf figure 21c. In contrast to Case A15c1, in which breakdown to small
scales is initiated in the inner layers adjacent to stable streaks, the main
mechanism of breakdown in Case C50c2 is the disintegration of the meandering
streak in the outer layer.

##  7 Conclusions and outlook

We have investigated the transition to turbulence in the bottom boundary layer
beneath a solitary wave by means of a simple parallel model taking into
account finite amplitude perturbations. The study consists of two steps
addressing the receptivity and breakdown stages of transition. In the
receptivity step, the most "dangerous" disturbances to which the boundary
layer shows the strongest response are found using a linear input-output
framework. In this framework, the perturbations are modelled as deterministic
body forces. The focus is in particular on early times prior to the flow
reversal. The optimal excitation per energy input was found to concentrate on
cross-stream components, which are arranged as streamwise-constant counter-
rotating rotational cells. These cells can be either steady or oscillate at
frequencies close to the effective wave frequency. This optimally-arranged
transverse forces introduce counter-rotating rollers that mix the streamwise
momentum of the flow and introduce energetic streamwise-constant streaks via
the lift-up effect. We have then selected a representative case with steady
streamwise-constant configuration at a spanwise wavenumber ($\beta=1.5$) to
seed small-amplitude rollers into nonlinear equations. As in the linear case,
the dynamics of the rollers are completely decoupled from the base flow and
the wave, hence they remain steady throughout the event. Optimally-arranged
steady rollers were found to amplify the energy of the streaks with a factor
proportional to $\Rey_{\delta}^{2}$. Increasing the amplitude of rollers
$\tilde{w}_{max}$ (58) leads to streaks with increased asymmetry, where low-
speed streaks become narrower and elevate into higher flow regions.

In the analysis of the breakdown step, we have first investigated the linear
secondary stability of perturbed boundary layers to identify the unstable
regions beneath the wave. To this end, we employed a quasi-static assumption,
which allows a separate stability analysis at each phase using the frozen base
flow. Two different streak instabilities were observed, which we denoted as
“inner” and “outer” instabilities after the location of their respective
critical layers, a naming convention suggested by Vaughan & Zaki (2011) for
flat-plate boundary layers. The inner instabilities have varicose symmetry and
are fed on the vertical shear, thus they have critical layers near to the
wall. They are activated in the APG stage at the same phases with the two-
dimensional instabilities of the baseline unperturbed flow. Compared to the
baseline instabilities, the inner instabilities have reduced growth rates due
to negative production driven by spanwise shear and enhanced dissipation in
two-dimensional mode shapes. The inner instabilities are therefore stabilizing
and can delay the transition to turbulence or completely suppress it. The
damping effect is strongest in streaks generated by rollers with magnitude
$\tilde{w}_{\max}\approx 2.8/\Rey_{\delta}$. In contrast to inner
instabilities, outer instabilities were found to be very unstable. They are of
sinuous nature and develop around the lifted low-speed streaks in the outer
region. These instabilities are driven by the spanwise shear of the base flow.
Therefore, they are only active when the low-speed streaks are significantly
elevated, which is achieved when the amplitude of the streaks $A_{s}$ (57)
exceeds 15% of the local free-stream velocity at the phase. This can occur
already in the FPG stage if the roller-perturbations are strong. Therefore,
outer instabilities can lead to a subcritical bypass transition at this stage.
The bifurcation point from inner to outer instabilities depends on the roller
magnitude and Reynolds number, and is found to be at $\tilde{w}_{\max}\approx
3.8/\Rey_{\delta}$.

In the final step of our analysis, the results of secondary stability analysis
were verified by means of DNS. We have observed a specific energy level above
which breakdown to turbulence occurred in all considered cases. Using this
empirical threshold, flow-state diagrams were generated. At a particular
phase, the state of the flow, i.e., laminar or turbulent, depends on Reynolds
number ($\Rey_{\delta}$), the roller amplitude ($\tilde{w}_{max}$) and the
initial amplitude of the tertiary perturbation in the secondary instability
mode. The state diagrams showed the damping effect of streaks more clearly,
e.g., the laminar zone protrudes deep into the turbulent zone for moderate-
amplitude perturbations. For instance, for the case
$\tilde{w}_{\max}=2.8/\Rey_{\delta}$, the damping mechanism can keep the flow
laminar up to very high Reynolds numbers such as $\Rey_{\delta}=4000$. These
observations suggest that the classification of flow states should at least
include an additional measure for environment perturbations. Previous
Reynolds-number based classifications are not satisfactory.

We have investigated the effect of finite amplitude perturbations on the
transition of a SWBL using an idealized deterministic model, which allows
generation of streaks in a controlled setting. A possible future direction is
extending the work to a more natural configuration, in which the ambient
turbulence and its penetration into boundary layer are considered. In this
model, streamwise vortices and streaks will evolve in a stochastic setting.
Depending on streak amplitudes, four possible transition scenarios are
anticipated: (i) orderly transition when streaks have negligible influence;
(ii) delayed transition under low- to moderate-amplitude ambient turbulence,
where inner instabilities on moderate-amplitude streaks dominate the APG
stage; (iii) bypass-transition under high-amplitude ambient turbulence, where
outer instabilities broke streaks into turbulent spots, which then grow, merge
and occupy the whole boundary layer; (iv) mixed transition, where any of the
prior transition mechanisms can occur at different parts of the boundary
layer. The mixed transition can occur in particular when the amalgamation
timescale of turbulent spots is slow. In this case, other transition
mechanisms can take place in laminar regions surrounding spots, e.g.,
turbulent spots and orderly spanwise rollers coexisted in the APG stage in
Sumer et al. (2010). Only after full assessment of the amalgamation timescale,
it will be clear under which circumstances a complete bypass transition can
take place in a SWBL.

## Acknowledgements

The research reported here has been supported by a grant from Ministry of
Education of Singapore to National University of Singapore. The computational
work for this article was fully performed on resources of the National
Supercomputing Centre, Singapore (https://www.nscc.sg). We thank Mengqi Zhang
for providing the starting point of our linear adjoint code and for useful
comments.

##  Appendix A Derivation of the optimal forcing

For our time-dependent problem, the adjoint approach can be utilized using the
formal Lagrange method (Corbett & Bottaro, 2001; Tröltzsch, 2010). First, the
inner products are defined as

$\langle\mbox{\boldmath$a$},\mbox{\boldmath$b$}\rangle_{\Omega}=\frac{1}{2}\int_{0}^{\infty}\mbox{\boldmath$($}\mbox{\boldmath$a$}^{*}\cdot\mbox{\boldmath$b$})\mathrm{d}z+c.c.;~{}~{}~{}\langle\mbox{\boldmath$a$},\mbox{\boldmath$b$}\rangle_{\overline{\Omega}}=\frac{1}{2}\int_{-\infty}^{T_{f}}\int_{0}^{\infty}\mbox{\boldmath$($}\mbox{\boldmath$a$}^{*}\cdot\mbox{\boldmath$b$})\mathrm{d}z\mathrm{d}t+c.c.,$
(76)

where asterisk denotes complex-conjugated fields and c.c. stands for the
complex conjugate of the previous terms in the expression. Subsequently, we
associate the following Lagrangian functional to the problem

$\displaystyle\mathcal{L}(\mbox{\boldmath$\hat{q}$},\mbox{\boldmath$\hat{q}$}^{+},\mbox{\boldmath$\hat{f}$},\sigma)$
$\displaystyle:=\frac{E(\mbox{\boldmath$\hat{q}$}(T_{f}))}{E(\mbox{\boldmath$\hat{q}$}(T_{i}))}+\langle\mbox{\boldmath$\hat{q}$}^{+},L(t)\mbox{\boldmath$\hat{q}$}-C\mbox{\boldmath$\hat{f}$}\mathrm{e}^{\mathrm{i}\omega_{f}t}\rangle_{\overline{\Omega}}+\sigma(\langle\mbox{\boldmath$\hat{f}$},\mbox{\boldmath$\hat{f}$}\rangle_{\Omega}-1),$
(77)

where $\mbox{\boldmath$\hat{q}$}^{+}$ is the Lagrange multiplier in the form
of adjoint perturbation fields to impose state constraints, and $\sigma$ is
the Lagrange multiplier to constrain the force to unity magnitude. In the
Lagrangian (77), we have employed an instant $T_{i}$ at which
$\hat{q}(T_{i}):=\hat{q}_{0}$ to remove $\hat{q}_{0}$ from the derivation and
simplify the process. The first order optimality conditions for Lagrangian
$\mathcal{L}$ dictates that variation of $\mathcal{L}$ with respect to
forward, adjoint and control variables vanish identically (e.g. Gunzburger
(2002)), i.e.,

$\frac{\partial\mathcal{L}}{\partial{\mbox{\boldmath$\hat{q}$}}}\delta\mbox{\boldmath$\hat{q}$}+\frac{\partial\mathcal{L}}{\partial{\mbox{\boldmath$\hat{q}^{+}$}}}\delta\mbox{\boldmath$\hat{q}^{+}$}+\frac{\partial\mathcal{L}}{\partial{\mbox{\boldmath$\hat{f}$}}}\delta\mbox{\boldmath$\hat{f}$}+\frac{\partial\mathcal{L}}{\partial\sigma}\delta\sigma=0,$
(78)

where the directional variation is defined as, e.g., for the arbitrary
variation $\delta\mbox{\boldmath$\hat{q}$}$ in state space,

$\frac{\partial\mathcal{L}}{\partial{\mbox{\boldmath$\hat{q}$}}}\delta\mbox{\boldmath$\hat{q}$}=\lim\limits_{\epsilon\rightarrow
0}\frac{\mathcal{L}(\mbox{\boldmath$\hat{q}$}+\epsilon\delta\mbox{\boldmath$\hat{q}$},\mbox{\boldmath$\hat{q}$}^{+},\mbox{\boldmath$\hat{f}$},\sigma)-\mathcal{L}(\mbox{\boldmath$\hat{q}$},\mbox{\boldmath$\hat{q}$}^{+},\mbox{\boldmath$\hat{f}$},\sigma)}{\epsilon}.$
(79)

Setting the variations
$\delta\mbox{\boldmath$\hat{q}$}^{+}=\delta\mbox{\boldmath$\hat{f}$}=0$, and
letting $\mbox{\boldmath$\hat{q}$}^{\prime}$ vary freely yields
$\mathcal{L}_{\mbox{\boldmath$\hat{q}$}}(\mbox{\boldmath$\hat{q}$}^{\prime})=0$.
These equations are manipulated by utilizing integration by parts in space and
time as many times as necessary until all differential operators on state
fields are moved on to adjoint fields. The resulting boundary integrals in
this process are eliminated by utilizing the homogeneous boundary conditions
of OSS equations.

Variation of the Lagrangian with respect to each component of the forcing
vector should vanish as a result of optimality condition in (78). Enforcing
this stationarity condition yields for the streamwise component

$\displaystyle
2\frac{\partial\mathcal{L}}{\partial\hat{f}_{u}}\delta\hat{f}_{u}$
$\displaystyle=2\sigma\int_{0}^{\infty}\hat{f}_{u}^{*}\delta\hat{f}_{u}\mathrm{d}z+\int_{-\infty}^{T_{f}}\int_{0}^{\infty}(\hat{w}^{+})^{*}\mathrm{i}\alpha\frac{\partial\delta\hat{f}_{u}}{\partial
z}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}z\mathrm{d}t$
$\displaystyle-\int_{-\infty}^{T_{f}}\int_{0}^{\infty}(\hat{\eta}^{+})^{*}\mathrm{i}\beta\delta\hat{f_{u}}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}z\mathrm{d}t+c.c.$
$\displaystyle=2\sigma\int_{0}^{\infty}\hat{f}_{u}^{*}\delta\hat{f}_{u}\mathrm{d}z+\left.\int_{-\infty}^{T_{f}}\left((\hat{w}^{+})^{*}\mathrm{i}\alpha\delta\hat{f}_{u}\right)\right|_{0}^{\infty}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}t$
$\displaystyle-\int_{-\infty}^{T_{f}}\int_{0}^{\infty}\frac{\partial(\hat{w}^{+})^{*}}{\partial
z}\mathrm{i}\alpha\delta\hat{f}_{u}\mathrm{d}z\mathrm{d}t-\int_{-\infty}^{T_{f}}\int_{0}^{\infty}(\hat{\eta}^{+})^{*}\mathrm{i}\beta\delta\hat{f_{u}}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}z\mathrm{d}t+c.c.$
$\displaystyle=\int_{0}^{\infty}\delta\hat{f_{u}}\left(2\sigma\hat{f}_{u}^{*}+\int_{-\infty}^{T_{f}}\left(-\mathrm{i}\alpha\frac{\partial(\hat{w}^{+})^{*}}{\partial
z}-\mathrm{i}\beta(\hat{\eta}^{+})^{*}\right)\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}t\right)\mathrm{d}z+c.c.=0,$

where we have made use of the Green’s identity and the homogeneous adjoint
boundary conditions $\hat{w}^{+}(0)=\hat{w}^{+}(z\rightarrow\infty)=0$. As the
variation $\delta\hat{f_{u}}$ is a free variable, the optimality condition
holds only if

$2\sigma\hat{f}_{u}^{*}+\int_{-\infty}^{T_{f}}\left(-\mathrm{i}\alpha\frac{\partial(\hat{w}^{+})^{*}}{\partial
z}-\mathrm{i}\beta(\hat{\eta}^{+})^{*}\right)\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}t=0.$

Manipulating the complex conjugates, we obtain the following expression for
the streamwise component of the optimal force

$\hat{f}_{u}^{opt}:=\hat{f}_{u}=-\frac{1}{2\sigma}\int_{-\infty}^{T_{f}}\left(\mathrm{i}\alpha\frac{\partial\hat{w}^{+}}{\partial
z}+\mathrm{i}\beta\hat{\eta}^{+}\right)\mathrm{e}^{-\mathrm{i}\omega_{f}t}\mathrm{d}t.$
(80)

The spanwise component of the optimal force is derived in a similar way:

$\displaystyle
2\frac{\partial\mathcal{L}}{\partial\hat{f}_{v}}\delta\hat{f}_{v}$
$\displaystyle=2\sigma\int_{0}^{\infty}\hat{f}_{v}^{*}\delta\hat{f}_{v}\mathrm{d}z+\int_{-\infty}^{T_{f}}\int_{0}^{\infty}(\hat{w}^{+})^{*}\mathrm{i}\beta\frac{\partial\delta\hat{f}_{v}}{\partial
z}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}z\mathrm{d}t$
$\displaystyle+\int_{-\infty}^{T_{f}}\int_{0}^{\infty}(\hat{\eta}^{+})^{*}\mathrm{i}\alpha\delta\hat{f_{v}}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}z\mathrm{d}t+c.c.$
$\displaystyle=2\sigma\int_{0}^{\infty}\hat{f}_{v}\delta\hat{f}_{v}\mathrm{d}z+\left.\int_{-\infty}^{T_{f}}(\hat{w}^{+})^{*}\mathrm{i}\beta\delta\hat{f}_{v}\right|_{0}^{z_{\max}}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}t$
$\displaystyle-\int_{-\infty}^{T_{f}}\int_{0}^{\infty}\frac{\partial(\hat{w}^{+})^{*}}{\partial
z}\mathrm{i}\beta\delta\hat{f}_{v}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}z\mathrm{d}t+\int_{-\infty}^{T_{f}}\int_{0}^{\infty}(\hat{\eta}^{+})^{*}\mathrm{i}\alpha\delta\hat{f_{v}}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}z\mathrm{d}t$
$\displaystyle=\int_{0}^{\infty}\delta\hat{f_{v}}\left(2\sigma\hat{f}_{v}^{*}+\int_{-\infty}^{T_{f}}\left(-\mathrm{i}\beta\frac{\partial(\hat{w}^{+})^{*}}{\partial
z}+\mathrm{i}\alpha(\hat{\eta}^{+})^{*}\right)\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}t\right)\mathrm{d}z=0,$

which yields for the optimal spanwise force

$\hat{f}_{v}^{opt}=\frac{1}{2\sigma}\int_{-\infty}^{T_{f}}\left(-\mathrm{i}\beta\frac{\partial\hat{w}^{+}}{\partial
z}+\mathrm{i}\alpha\hat{\eta}^{+}\right)\mathrm{e}^{-\mathrm{i}\omega_{f}t}\mathrm{d}t.$
(81)

Furthermore, the vertical component is derived as follows

$\displaystyle
2\frac{\partial\mathcal{L}}{\partial\hat{f}_{w}}\delta\hat{f}_{w}$
$\displaystyle=2\sigma\int_{0}^{\infty}\hat{f}_{w}^{*}\delta\hat{f}_{w}\mathrm{d}z+\int_{-\infty}^{T_{f}}\int_{0}^{\infty}(\hat{w}^{+})^{*}k^{2}\delta\hat{f}_{w}\mathrm{d}z\mathrm{d}t$
$\displaystyle=\int_{0}^{\infty}\delta\hat{f_{w}}\left(2\sigma\hat{f}_{w}^{*}+\int_{-\infty}^{T_{f}}k^{2}(\hat{w}^{+})^{*}\mathrm{e}^{\mathrm{i}\omega_{f}t}\mathrm{d}t\right)\mathrm{d}z=0.$

Consequently, we obtain

$\hat{f}_{w}^{opt}=-\frac{1}{2\sigma}\int_{-\infty}^{T_{f}}k^{2}\hat{w}^{+}\mathrm{e}^{-\mathrm{i}\omega_{f}t}\mathrm{d}t.$
(82)

Finally, the variation with respect to $\sigma$ at a stationary point is

$\frac{\partial\mathcal{L}}{\partial\sigma}\delta\sigma=\delta\sigma(\langle\mbox{\boldmath$\hat{f}$},\mbox{\boldmath$\hat{f}$}\rangle_{\Omega}-1)=0,$

which restores the constraint equation for the amplitude of the forcing

$\langle\mbox{\boldmath$\hat{f}$},\mbox{\boldmath$\hat{f}$}\rangle_{\Omega}=1$
(83)

Equations (80)–(83) represent a closed system of equations for the three
forcing components and $\sigma$.

## References

  * Andersson et al. (1999) Andersson, P., Berggren, M. & Henningson, D. S. 1999 Optimal disturbances and bypass transition in boundary layers. Physics of Fluids 11 (1), 134–150.
  * Andersson et al. (2001) Andersson, P., Brandt, L., Bottaro, A. & Henningson, D. S. 2001 On the breakdown of boundary layer streaks. Journal of Fluid Mechanics 428, 29–60.
  * Barkley et al. (2008) Barkley, D., Blackburn, H. M. & Sherwin, S. J. 2008 Direct optimal growth analysis for timesteppers. International Journal for Numerical Methods in Fluids 57 (9), 1435–1458.
  * Blondeaux et al. (2012) Blondeaux, P., Pralits, J. & Vittori, G. 2012 Transition to turbulence at the bottom of a solitary wave. Journal of Fluid Mechanics 709, 396–407.
  * Bolis (2013) Bolis, A. 2013 Fourier spectral/hp element method: investigation of time-stepping and parallelisation strategies. PhD thesis, Imperial College London.
  * Boyd (2001) Boyd, J. P. 2001 Chebyshev and Fourier spectral methods. Courier Corporation.
  * Brandt (2014) Brandt, L. 2014 The lift-up effect: the linear mechanism behind transition and turbulence in shear flows. European Journal of Mechanics-B/Fluids 47, 80–96.
  * Butler & Farrell (1992) Butler, K. M. & Farrell, B. F. 1992 Three-dimensional optimal perturbations in viscous shear flow. Physics of Fluids 4 (8), 1637–1650.
  * Cantwell et al. (2015) Cantwell, C. D., Moxey, D., Comerford, A., Bolis, A., Rocco, G., Mengaldo, G., De Grazia, D., Yakovlev, S., Lombard, J. E., Ekelschot, D., Jordi, B., Xu, H., Mohamied, Y., Eskilsson, C., Nelson, B., Vos, P., Biotto, C., Kirby, R. M. & Sherwin, S. J. 2015 Nektar++: An open-source spectral/hp element framework. Computer Physics Communications 192, 205–219.
  * Carstensen et al. (2010) Carstensen, S., Sumer, B. M. & Fredsøe, J. 2010 Coherent structures in wave boundary layers. Part 1. Oscillatory motion. Journal of Fluid Mechanics 646, 169–38.
  * Chen & Kirchner (1971) Chen, C. F. & Kirchner, R. P. 1971 Stability of time-dependent rotational couette flow part 2. stability analysis. Journal of Fluid Mechanics 48 (2), 365–384.
  * Corbett & Bottaro (2001) Corbett, P. & Bottaro, A. 2001 Optimal linear growth in swept boundary layers. Journal of Fluid Mechanics 435 (0), 1–23.
  * Cossu & Brandt (2002) Cossu, C. & Brandt, L. 2002 Stabilization of tollmien–schlichting waves by finite amplitude optimal streaks in the blasius boundary layer. Physics of Fluids 14 (8), L57–L60.
  * Cossu & Brandt (2004) Cossu, C. & Brandt, L. 2004 On tollmien–schlichting-like waves in streaky boundary layers. European Journal of Mechanics-B/Fluids 23 (6), 815–833.
  * Durbin (2017) Durbin, P. A. 2017 Perspectives on the Phenomenology and Modeling of Boundary Layer Transition. Flow, Turbulence and Combustion 99 (1), 1–23.
  * Grimshaw (1970) Grimshaw, R 1970 The solitary wave in water of variable depth. Journal of Fluid Mechanics 42 (3), 639–656.
  * Gunzburger (2002) Gunzburger, M. 2002 Perspectives in Flow Control and Optimization. Society for Industrial and Applied Mathematics.
  * Gustavsson (1991) Gustavsson, L. H. 1991 Energy growth of three-dimensional disturbances in plane Poiseuille flow. Journal of Fluid Mechanics 224, 241.
  * Hack & Zaki (2014) Hack, M. J. P. & Zaki, T. A. 2014 Streak instabilities in boundary layers beneath free-stream turbulence. Journal of Fluid Mechanics 741, 280–315.
  * Hernon et al. (2007) Hernon, D., Walsh, E. J. & McELIGOT, D. M. 2007 Experimental investigation into the routes to bypass transition and the shear-sheltering phenomenon. Journal of Fluid Mechanics 591, 21\.
  * Hunt et al. (1988) Hunt, J. C. R., Wray, A. A. & Moin, P. 1988 Eddies, streams, and convergence zones in turbulent flows. In Proceedings of the 1988 Summer Program, pp. 193–208.
  * Jacobs & Durbin (2001) Jacobs, R. G. & Durbin, P. A. 2001 Simulations of bypass transition. Journal of Fluid Mechanics 428 (0), 185–212.
  * Jovanović & Bamieh (2005) Jovanović, M. R. & Bamieh, B. 2005 Componentwise energy amplification in channel flows. Journal of Fluid Mechanics 534, 145–183.
  * Karniadakis (1990) Karniadakis, G. E. 1990 Spectral element-Fourier methods for incompressible turbulent flows. Computer Methods in Applied Mechanics and Engineering 80 (1), 367–380.
  * Karniadakis et al. (1991) Karniadakis, G. E., Israeli, M. & Orszag, S. A. 1991 High-order splitting methods for the incompressible navier-stokes equations. Journal of Computational Physics 97 (2), 414–443.
  * Karniadakis & Sherwin (2005) Karniadakis, G. E. & Sherwin, S. J. 2005 Spectral/hp methods for computational fluid dynamics. Oxford University Press.
  * Kirby & Sherwin (2006) Kirby, R. M. & Sherwin, S. J. 2006 Aliasing errors due to quadratic nonlinearities on triangular spectral /hp element discretisations. Journal of Engineering Mathematics 56 (3), 273–288.
  * Klebanoff et al. (1962) Klebanoff, P. S., Tidstrom, K. D. & Sargent, L. M. 1962 The three-dimensional nature of boundary-layer instability. Journal of Fluid Mechanics 12 (1), 1–34.
  * Klettner & Eames (2012) Klettner, C. A. & Eames, I. 2012 The laminar free surface boundary layer of a solitary wave. Journal of Fluid Mechanics 696, 423–433.
  * Landahl (1980) Landahl, M. T. 1980 A note on an algebraic instability of inviscid parallel shear flows. Journal of Fluid Mechanics 98, 243–251.
  * Liu & Orfila (2004) Liu, P. L. F. & Orfila, A. 2004 Viscous effects on transient long-wave propagation. Journal of Fluid Mechanics 520, 83–92.
  * Liu et al. (2007) Liu, P. L. F., Park, Y. S. & Cowen, E. A. 2007 Boundary layer flow and bed shear stress under a solitary wave. Journal of Fluid Mechanics 574, 449.
  * Liu et al. (2008a) Liu, Y., Zaki, T. A. & Durbin, P. A. 2008a Boundary-layer transition by interaction of discrete and continuous modes. Journal of Fluid Mechanics 604, 199–233.
  * Liu et al. (2008b) Liu, Y., Zaki, T. A. & Durbin, P. A. 2008b Floquet analysis of secondary instability of boundary layers distorted by klebanoff streaks and tollmien–schlichting waves. Physics of Fluids 20 (12), 124102\.
  * Luchini & Bottaro (2014) Luchini, P. & Bottaro, A. 2014 Adjoint Equations in Stability Analysis. Annual Review of Fluid Mechanics 46 (1), 493–517.
  * Mans et al. (2007) Mans, J., de Lange, H. C. & van Steenhoven, A. A. 2007 Sinuous breakdown in a flat plate boundary layer exposed to free-stream turbulence. Physics of Fluids 19 (8), 088101–5.
  * Matsubara & Alfredsson (2001) Matsubara, M. & Alfredsson, P. H. 2001 Disturbance growth in boundary layers subjected to free-stream turbulence. Journal of Fluid Mechanics 430 (0), 149–168.
  * Morkovin (1969) Morkovin, M. V. 1969 On the many faces of transition. In Viscous drag reduction, pp. 1–31. Springer.
  * Munk (1949) Munk, W. H. 1949 The solitary wave theory and its application to surf problems. Annals of the New York Academy of Sciences 51 (3), 376–424.
  * Narasimha (1985) Narasimha, R. 1985 The laminar-turbulent transition zone in the boundary layer. Progress in Aerospace Sciences 22 (1), 29–80.
  * Ozdemir et al. (2013) Ozdemir, C. E., Hsu, T.-J. & Balachandar, S. 2013 Direct numerical simulations of instability and boundary layer turbulence under a solitary wave. Journal of Fluid Mechanics 731, 545–578.
  * Phillips (1969) Phillips, O. M. 1969 Shear-flow turbulence. Annual Review of Fluid Mechanics 1 (1), 245–264.
  * Reddy et al. (1998) Reddy, S. C., Schmid, P. J., Baggett, J. S. & Henningson, D. S. 1998 On stability of streamwise streaks and transition thresholds in plane channel flows. Journal of Fluid Mechanics 365 (0), 269–303.
  * Ricco et al. (2011) Ricco, P., Luo, J. & Wu, X. 2011 Evolution and instability of unsteady nonlinear streaks generated by free-stream vortical disturbances. Journal of Fluid Mechanics 677, 1–38.
  * Rocco (2014) Rocco, G. 2014 Advanced instability methods using spectral/hp discretisations and their applications to complex geometries. PhD thesis, Imperial College London.
  * Sadek (2015) Sadek, M. 2015 Mechanisms for transition to turbulence in the bottom boundary layer under a surface solitary wave. PhD thesis, Cornell University.
  * Sadek et al. (2015) Sadek, M. M., Parras, L., Diamessis, P. J. & LIU, P. L. F. 2015 Two-dimensional instability of the bottom boundary layer under a solitary wave. Physics of Fluids 27 (4), 044101–25.
  * Schmid & Brandt (2014) Schmid, P. J. & Brandt, L. 2014 Analysis of Fluid Systems: Stability, Receptivity, Sensitivity. Applied Mechanics Reviews 66 (2), 021003–21.
  * Schmid & Henningson (2001) Schmid, P. J. & Henningson, D. S. 2001 Stability and Transition in Shear Flows, Applied Mathematical Sciences, vol. 142. New York, NY: Springer New York.
  * Sumer et al. (2010) Sumer, B. M., Jensen, P. M., Sørensen, L. B., Fredsøe, J., Liu, P. L. F. & Carstensen, S. 2010 Coherent structures in wave boundary layers. Part 2. Solitary motion. Journal of Fluid Mechanics 646, 207–25.
  * Tanaka et al. (2012) Tanaka, Hitoshi, Winarta, Bambang, Suntoyo & Yamaji, Hiroto 2012 Validation of a new generation system for bottom boundary layer beneath solitary wave. Coastal Engineering 59 (1), 46–56.
  * Trefethen et al. (1993) Trefethen, L. N., Trefethen, A. E., Reddy, S. C. & Driscoll, T. A. 1993 Hydrodynamic stability without eigenvalues. Science 261 (5121), 578–584.
  * Tröltzsch (2010) Tröltzsch, F. 2010 Optimal control of partial differential equations: theory, methods, and applications. American Mathematical Society.
  * Tuckerman & Barkley (2000) Tuckerman, L. S. & Barkley, D. 2000 Bifurcation analysis for timesteppers. In Numerical methods for bifurcation problems and large-scale dynamical systems, pp. 453–466. Springer.
  * Vaughan & Zaki (2011) Vaughan, N. J. & Zaki, T. A. 2011 Stability of zero-pressure-gradient boundary layer distorted by unsteady Klebanoff streaks. Journal of Fluid Mechanics 681, 116–153.
  * Verschaeve & Pedersen (2014) Verschaeve, J. C. G. & Pedersen, G. K. 2014 Linear stability of boundary layers under solitary waves. Journal of Fluid Mechanics 761, 62–104.
  * Verschaeve et al. (2017) Verschaeve, J. C. G., Pedersen, G. K. & Tropea, C. 2017 Non-modal stability analysis of the boundary layer under solitary waves. Journal of Fluid Mechanics 836, 740–772.
  * Vittori & Blondeaux (2008) Vittori, G. & Blondeaux, P. 2008 Turbulent boundary layer under a solitary wave. Journal of Fluid Mechanics 615, 433\.
  * Vittori & Blondeaux (2011) Vittori, G. & Blondeaux, P. 2011 Characteristics of the boundary layer at the bottom of a solitary wave. Coastal Engineering 58 (2), 206–213.
  * Von Kerczek & Davis (1974) Von Kerczek, C. & Davis, S. H. 1974 Linear stability theory of oscillatory stokes layers. Journal of Fluid Mechanics 62 (4), 753–773.
  * Vos et al. (2011) Vos, P. E. J., Eskilsson, C., Bolis, A., Chun, S., Kirby, R. M. & Sherwin, S. J. 2011 A generic framework for time-stepping partial differential equations (PDEs): general linear methods, object-oriented implementation and application to fluid problems. International Journal of Computational Fluid Dynamics 25 (3), 107–125.
  * Waleffe (1995) Waleffe, F. 1995 Hydrodynamic stability and turbulence: Beyond transients to a self-sustaining process. Studies in Applied Mathematics 95 (3), 319–343.
  * Weideman & Reddy (2000) Weideman, J. A. C. & Reddy, S. C. 2000 A MATLAB differentiation matrix suite. ACM Transactions on Mathematical Software 26 (4), 465–519.