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

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.


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  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.
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 nearwall 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 δ 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 twodimensional 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 freeshear layers. In contrast, in high-noise cases, where the streaks should be strongest, vortex rollers broke into Λ-shaped vortices, hence the transition was reminiscent of Ktype 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 Re δ = 1000. These sporadic features spread to earlier phases with increasing Reynolds number. At the highest Reynolds number achieved (Re δ = 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 , 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 streamwiseconstant 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  or modify the secondary modes of transition when they have modest energy (Ozdemir et al. 2013). Furthermore, we anticipate that modestamplitude 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.

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) where c * w = √ 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 with the maximum velocity 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 is the effective wave frequency. Using the wave frequency and kinematic viscosity the Stokes length is defined as the boundary-layer scale of the problem. Equation (2.4) neglects the first and higherorder terms in H * /h * . Therefore, the model is relevant only for H * /h * → 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 * → 0 and δ * /h * → 0 provides the advantage of reducing the parameter space of the problem to one, the Reynolds number: 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, The momentum equation for the irrotational streamwise velocity at the bottom can be expressed in local temporal frame approximation as 2 Re δ The free-stream pressure gradient is calculated using (2.4) and (2.9) and reads (2.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 where U = u r +u 0 is the total laminar velocity containing both rotational and irrotational components. Equation (2.11) is supplemented with the boundary conditions U (z = 0, t) = 0 and U (z → ∞, t) = u 0 , and we specifiy the initial condition U (z, −∞) = 0. The solution of (2.11) is shown in figure 2c. This approximate model of SWBL is theoretically solved (Liu & Orfila 2004) and adapted in experimental Tanaka et al. 2012) and numerical (Ozdemir et al. 2013;Sadek et al. 2015;Verschaeve et al. 2017) studies.

Optimal disturbances and the flow response
The time-dependent streamwise velocity U (z, t) in (2.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 where α and β are streamwise and spanwise wavenumbers, respectively, and ω 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.

Methodology
We apply a Fourier ansatz in the homogeneous x and y directions for the perturbation velocity and pressure [ũ,ṽ,w,p](x, y, z, t) = Re{ [û,v,ŵ,p](z, t)e i(αx+βy) }. (3.2) In the linear regime, each Fourier mode is excited by the corresponding harmonic force at the same spatial wavenumber 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 −π, the base flow is vanishingly small, and the forcing brings the stagnant flow to a periodic state, which is given by the set whereŵ andη are the vertical velocity and vertical vorticity modes associated with the wavenumber pair (α, β),∆ = k 2 − ∂ 2 /∂z 2 represents the semi-discretised Laplacian operator, k 2 = α 2 + β 2 , andĝ w andĝ η are the external driving forces containing the When the wave arrives, the flow is no longer periodic, and the perturbation equations during the wave event becomes (3.14) The following compact notation is used for periodic and temporal OSS system of equations The response of the flow to an excitation at a wavenumber pair (α, β) is measured by the perturbation kinetic energy We can expressû andv in terms ofŵ andη, and obtain (Schmid & Henningson 2001) E(q) := 1 2k 2 ∞ 0 (k 2 |ŵ| 2 + | ∂ŵ ∂z | 2 + |η| 2 )dz. (3.17) We look for the most dangerous perturbations initiating the strongest response in the linear SWBL. This is equivalent to finding an optimal controlf opt (z; α, β, ω f , T f , Re δ ), which yields the maximum energy amplification at a terminal time t = T f per initial energy input E(q o ). This is found by solving the constrained optimization problem where G f is the largest response or gain. The optimization problem (3.18) 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 (3.18) 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 configurationf opt satisfying L oqo = Cf opt , The reader is directed to Schmid & Henningson (2001) for a thorough derivation of equations (3.21) and (3.22). Furthermore, (3.20) corresponds to the expressions to calculate the optimal forcing configuration using the adjoint fields, i.e., where σ is a Lagrange multiplier, cf. Appendix A for the details of the derivation of (3.26)-(3.28). The optimization problem in (3.18) is now transformed to a set of equations with (3.4)-(3.7) and (3.10)-(3.14) being the state or forward equations, (3.21)-(3.25) being the adjoint equations, and (3.26)-(3.28) 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 off opt and iterates over the following successive steps: (i) calculation ofq o using (3.4)-(3.7); (ii) forward-in-time integration of the state equations (3.10)-(3.14); (ii) backward-in-time integration of the adjoint equations in (3.21) and (3.25); (iii) updating the control terms with the available adjoint fields using (3.26)-(3.28) and 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 ξ ∈ [−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 ∈ [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π. At this phase the effect of the wave is negligible. The Crank-Nicolson scheme is employed for time integration. The time step size is δt = π/480. A sensitivity analysis confirmed that the selected spatial and temporal resolutions are sufficient.

Linear response of the flow
In this section, we study the linear response of the flow to the optimal perturbations at different α, β, ω f , T f and Re δ . Figure 3 shows the maximum gain among all excitation frequencies for each wavenumber pair (α, β) at a moderate Reynolds number Re δ = 2000. Re δ = 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 streamwiseconstant (α = 0, β = 0) excitation in the FPG stage and has a very weak response to two-dimensional α = 0, β = 0 and oblique α = 0, β = 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 4 further demonstrates the response of the flow to streamwise-constant excitation on a β − ω f plane at several terminal times (T f = −π/3, −π/6, 0, and π/6) at Re δ = 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 β ∈ [1.5, 2.5] and ω f ∈ [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 = π/6), cf. figure 4d.
In order to analyze the scaling of the governing equations with Reynolds number in the case of streamwise-constant excitation, we introduce the transformations and substitude to streamwise-constant OSS equations (3.10) and (3.11), where the terms with α vanish, i.e., In the scaled streamwise-constant setting, the velocity components becomes Therefore, the evolution of cross-stream momentum by the velocity components v and w is embedded into (3.30) and the evolution of streamwise momentum by u is linked to (3.31). As shown in (3.30) 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.,v =v o e iω f t ,ŵ =ŵ o e iω f t .
(3.33) Therefore, the increase in G f is solely due to intrinsic amplification ofû by the boundary layer. Equation (3.30) suggests that introducing w rendered the cross-stream momentum balance independent of Reynold number, while (3.31) indicates that the streamwise forcing is one order lower (O(1/Re δ )) 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., (3.34) Consequently, the streamwise forcing in (3.31) can be neglected for Re δ 1, and the evolution of u becomes independent of Reynolds number. Figure 5 validates these Reynolds-number scalings using the numerical results for the case T f = 0, α = 0, β = 1.5, and ω 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 2 δ and the terminal vertical velocity scaled with Re δ collapse for different Reynolds numbers, cf. figures 5b and 5c.
We now turn to input and output configurations. The optimal steady streamwiseconstant forcing configuration f opt (α = 0, β = 1.5, ω 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 It is observed that the steady forcing is organized as counter-rotating cells (0,f opt v ,f opt w ), which induce steady rollers (0,ṽ o ,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ũ and positivew 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 (3.30) and (3.31), η and w are of the same order in Reynolds number. So, a transverse steady forcing with an amplitude of O(1/Re 2 δ ) will induce steady rollers of strength O(1/Re δ ). 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 streakroller system, where O(1) streaks synthetically forced by the steady rollers of magnitude O(1/Re δ ). 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).

Nonlinear streaks
When the perturbations reach appreciable amplitudes, nonlinear effects should be taken into account. We showed in § 3.2 that the cases with β = 0, α = 0 and ω 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 crossstream 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 where A o is a small forcing magnitude with A o 1 and ∆ 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., where the small streamwise forcing A o f opt u is neglected. As the evolution of the crossstream 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 and define the transformations Introducing the transformed variables to the roller equations (4.9) and to the streak equation Transforming the nonlinear governing equations from (4.1)-(4.4) to (4.7)-(4.10) reduces the parameter space of the problem from two, Re δ and A 0 , to one, A, which can be considered now as the effective amplitude of the excitation. We reiterate that this oneparameter model is only applicable in the range of Re δ 1, where optimal forcing configuration f opt does not depend on Re δ 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 ∈ [0, L z = 20] in the vertical direction, and to y ∈ [0, Ly = 2π/β] 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).
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 β 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 β = 1.5 shows good performance for time horizons corresponding to the strongest amplifications (T f = 0, π/6). Therefore, we will merely consider roller perturbations induced by optimal forcing f opt (α = 0, β = 1.5, ω 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 (4.12) Figure 7a shows the variation of the roller magnitudes with respect to effective forcing amplitude A. The amplitudes are presented in a Re δ -independent scaling, i.e.,w max Re δ = A ≈ 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 magnitudesw max = 2.8/Re δ , 9.51/Re δ and 18.78/Re δ . 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 c s = 0.152) and by Andersson et al. (2001) (A c s = 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.
The streaky fields induced by linearly optimal steady streamwise-constant forcing are  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 ≈ ψ 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 are the integrated kinetic energy at a time t and the initial energy, respectively. Although the streak-roller system given by [ũ, 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ṽ 2 o andw 2 o are two orders lower in Reynolds number thanũ 2 and therefore, can be neglected in sufficiently high Reynolds numbers, i.e., E V,u ≈ E V , where E V,u is the integrated streamwise kinetic energy. Consequently, if we define a gain with similarity variables then the quadratic dependency of G f on Reynolds number can be shown explicitly (4.18) Figure 9 shows the gains for A = 15, 50 and 100 at Re δ = 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 ≈ 0.9, and then decay with the flow reversal.

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 streamwiseinvariant, 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 (4.13). We consider three-dimensional tertiary perturbations of the form [u , v , w , p ](x, y, z, t) = Re{[û ,v ,ŵ ,p ](y, z, t)e i(αx− t 0 ω(τ )dτ ) }, where α are real wavenumbers, and ω = ω r + iω 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 These equations are complemented with boundary conditionsû =v =ŵ = 0 on z = 0 and z → ∞. The system of equations (5.2)-(5.5) can be written shortly as L (U (t))q = ω(t)q , whereq (y, z, t) = [û ,v ,ŵ ,p ](y, z, t) . Using the quasi-static approximation, an eigenvalue problem is defined by freezing the operator L at a temporal station t = t s and solving for ω(t s ) andq (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.
As in § 4, we only consider nonlinear streaks induced by streamwise-constant timeinvariant excitation f opt (α = 0, β = 1.5, ω f = 0, T f = 0, Re δ ). Figure 10b shows the maximum leading eigenvalues along all streamwise wavenumbers, ω max for varying excitation amplitudes A, (4.5). The stability boundaries are demonstrated with thick contourlines in figure 10b. A slightly positive value ω max i /Re δ = 0.0001 is employed, as the exact neutral curve (ω max i = 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 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 Re δ = 2000 and Re δ = 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 ≈ 1), the growth rate in A = 15

Eα(π)
outer instabilities inner instabilities Figure 11. The overall growth measured by the modal perturbation kinetic energy density Eα at t = π, cf. (5.6). The initial energy density is set to E0 = 1. The most unstable streamwise wavenumber (α) 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.
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, ω max i /Re δ , 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).
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 (5.1), the energy at a mode is expressed by where E 0 is the initial energy density at the mode. Figure 11 demonstrates the variation of E α at t = π with respect to A. In this figure, the most unstable α 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 <w max < 3.76/Re δ using figure 7a), and (iii) outer-instability regime (A > 20, orw max > 3.76/Re δ ). 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 (∂U s /∂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 ≈ 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 12 demonstrates the variation of growth rates with respect to streamwise wavenumbers α for the cases A = 0, 15, 50 and 100 calculated at Re δ = 2000. It is observed that the primary instabilities have relatively longer wavelength peaking at wavenumbers around α ≈ 0.45 (figure 12a), where the outer instabilities concentrate in shorter wavelengths, e.g., α ≈ β/2 ≈ 0.75 for A = 50 (figure 12c) and α ≈ β ≈ 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 α > 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 ≈ 1 and α ≈ 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 α ≈ 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 α ) at α = 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, α = 0.45); (A = 15, α = 0.35); (A = 50, α = 0.75); (A = 100, α = 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 c s : U s = c r = ω r /α, form. The eigenmodes concentrate in the critical layers, where they convect with local mean velocity U c s , cf. dashed lines in figure 13. Streaky boundary layers consist of spanwise and vertical shear layers associated with ∂U s /∂y and ∂U s /∂z, respectively. These are shown in figure 14 for two representative cases (A = 15, α = 0.35) and (A = 50, α = 0.75) at time t = π/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 outerinstability 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 where e V is the total perturbation kinetic energy, p V,y is the total production rate due to spanwise shear, p V,w is the total production rate due to vertical shear and d V is the total dissipation rate. Here we neglect the contributions related to base spanwise (V s ) and vertical (W s ) velocities asṽ o andw o are an order of magnitude smaller in Re δ . If we consider the perturbations associated with a single instability mode, we can express individual terms by λ y = 2π/β and E = 2 Re δ (û * û +v * v +ŵ * ŵ ),D = 2 Re δ (ε * ε +ζ * ζ +η * η ), (5.9) In the dissipation termD, the components of the perturbation vorticity [ε , ζ , η ] are employed. Now the growth rate of the instability can be expressed as  Figure 15 shows the production and dissipation fields for the cases A = 15, α = 0.35 and A = 50, α = 0.75 at time t = π/6. The total integrated values,p V,y ,p V,y ,d V , are also presented in the respective panels of the fields. The energy of the perturbations is normalized to unity, i.e., e 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. (5.10). In case A = 15, the shear-damping effect is stronger with |p V,y | being about 20% of p 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 lowamplitude streaks, Cossu & Brandt (2004) has reported inner TS-like instabilities with reduced growth rates. They suggested the negative contributions fromp V,y as the primary mechanism behind the stabilizing effect of low-amplitude streaks. This term vanishes in the unperturbed boundary layer (∂U/∂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 (ζ ), which yields lower dissipation rates, cf.D in (5.9). 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.  Andersson et al. (2001) for a outer streak instability in a ZPG boundary layer. Dotted line shows cr/u0(t) = 0.31, which is the phase velocity calculated by Cossu & Brandt (2004) for an inner (modified TS-like) instability. 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)).

Direct numerical simulations
In § 4, we analyzed the dynamics of streamwise-constant streaks in a two-dimensional domain (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 (α = 0, β = 1.5, ω f = 0, T f = 0, Re δ ).
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 = −π. 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 ] × [0, L y ] × [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).
The energy density in each streamwise mode α is calculated by Fourier transforming the velocity fields in the streamwise direction and integrating the respective energy in the Fourier modeû(α, y, z, t) over the domain and then normalizing it, i.e., Ly 0û * (α, y, z, t) ·û(α, y, z, t)dydz. (6.1) 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 (4.15), 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 α = 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 α = 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 α = 0.75 as in the Case A50c2, cf. figure 18c. There is also good agreement between DNS and LST in the cases containing innerinstabilities (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 nonnegligible 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 ≈ 1.3.
If we assume that LST results are always applicable and all instabilities are of inviscid nature with constant ω i /Re δ , then we can utilize the empirical threshold E c and the growth rates ω i /Re δ calculated at a specific Reynolds number (e.g. Re δ = 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,w max , Re δ , and the initial perturbation energy in the instability mode, E 0 . Figures 19 show the flow states with respect tow max and Re δ at t = 2/9π and t = π/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 tow max = 3.76/Re δ (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π. 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, 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 < Re δ < 1500 in figures 19a,b.
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 rangew max ≈ 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 ofw max = 2.8/Re δ , 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 19. State of the SWBL with respect to Reynolds number (Re δ ) and amplitude of steady-roller perturbations (wmax) are shown at two representative times in the APG stage for two different initial tertiary perturbations (E0). The rollers are induced by steady streamwise-constant forcing f opt (α = 0, β = 1.5, ω f = 0, T f = 0). (a) t = 2/9π, E0 = 10 −32 ; (b) t = 2/9π, E0 = 10 −20 ; (c) t = π/2, E0 = 10 −32 ; (d) t = π/2, E0 = 10 −20 . The red dashed lines demonstrate the boundary (A = 20) between inner and outer instabilities, which corresponds towmax = 3.76/Re δ (figure 7a). Figure 20 shows the breakdown of inner instability in Case A15c1. At t = 24/90π in figure 20a, we see at the center a low-speed streak making undulations in the downstream direction with wavenumber α = 0.35. Since the inner instability is of varicose nature, the undulations are symmetric with respect to the streak. Vortical structures around the lowspeed streak are also shown in the figure using a positive isosurface of Q-criterion (Hunt et al. 1988). Among several vortical features, Λ-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 Λ vortices developing on streak-modulated instability waves in ZPG boundary layers (Liu et al. 2008a). Later at t = 25/90π, 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π, cf. figure 20c.
The transition to turbulence in Case C50c2 is demonstrated in figure 21. Initially at t = 50/180π, 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 α = 0.75. Subsequently, at t = 52/180π, 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π, 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.

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 counterrotating 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 (β = 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 Re 2 δ . Increasing the amplitude of rollersw max (4.12) leads to streaks with increased asymmetry, where lowspeed 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 magnitudẽ w max ≈ 2.8/Re δ . 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 (4.11) exceeds 15% of the local free-stream velocity at the phase. This can occur already in the FPG stage if the rollerperturbations 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 atw max ≈ 3.8/Re δ .
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 (Re δ ), the roller amplitude (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 casew max = 2.8/Re δ , the damping mechanism can keep the flow laminar up to very high Reynolds numbers such as Re δ = 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. and control variables vanish identically (e.g. Gunzburger (2002) where the directional variation is defined as, e.g., for the arbitrary variation δq in state space, ∂L ∂q δq = lim →0 L(q + δq,q + ,f , σ) − L(q,q + ,f , σ) .

(A 4)
Setting the variations δq + = δf = 0, and lettingq vary freely yields Lq(q ) = 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 (A 3). Enforcing this stationarity condition yields for the streamwise component where we have made use of the Green's identity and the homogeneous adjoint boundary conditionsŵ + (0) =ŵ + (z → ∞) = 0. As the variation δf u is a free variable, the optimality condition holds only if Manipulating the complex conjugates, we obtain the following expression for the streamwise component of the optimal forcê f opt u :=f u = − 1 2σ