The objective of this paper is to establish a detailed map for the temporal instability of the bottom boundary layer (BBL) flow driven by a soliton-like wave induced pressure gradient in a U-tube-shaped tunnel, which serves as an approximation to the BBL under surface solitary waves. Both linear stability analysis and fully nonlinear two-dimensional simulations using high-order numerical methods have been carried out. The process of delineation of the stability regions as a function of boundary layer thickness-based Reynolds number of the temporally evolving base flow, *Re*_{δ}, consists of two parts. In the first part, we assess the lower limit of the *Re*_{δ} range within which the standard, quasi-steady, linear stability analysis is applicable when considering individual base flow profiles sampled during its transient evolution. Below this limit, transient linear stability analysis serves as a more accurate predictor of the stability properties of the base flow. In the second step, above the *Re*_{δ} limit where the BBL is determined to be linearly unstable, the base flow is further classified as unconditionally stable, conditionally unstable, or unconditionally unstable in terms of its sensitivity to the amplitude and the insertion time of perturbations. Two distinct modes of instability exist in this case: post- and pre-flow-reversal modes. At a moderate value of *Re*_{δ}, both modes are first observed in the wave deceleration phase. The post flow reversal mode dominates for relatively low *Re* and it is the one observed in Sumer *et al.* [“Coherent structures in wave boundary layers. Part 2. Solitary motion,” J. Fluid Mech. **646**, 207 (2010)]. For *Re* above a threshold value of the base flow in the unconditionally unstable regime, the pre-flow reversal mode, which is longer in wavelength than its post-reversal counterpart, becomes dominant. In the same regime, the threshold *Re*_{δ} value above which instability is observed in the acceleration phase of the wave is also identified. In this case, the base flow velocity profiles lack any inflection point, suggesting that the origin of such an instability is viscous. Finally, the lower *Re*_{δ} limit above which quasi-steady linear stability analysis is valid may be independently obtained by adapting to the surface solitary wave BBL, an instability criterion which links the average growth rate and wave event timescale, as previously proposed in studies of the instability of the interior of progressive and solitary internal waves.

## I. INTRODUCTION

Tsunami waves, which can be characterized as transient long waves, have a massively destructive impact on inundated coastal areas (e.g., 2004 Indian Ocean tsunami and 2011 Japan Tohoku tsunami). Two of the phenomena related to such events are the associated sediment deposition and sea bed erosion.^{1,2} Both of these phenomena bear crucial implications for the design of the off- and near-shore structures and for shoreline protection. These sediment-specific processes are mainly driven by the near-bed flow characteristics, i.e., turbulence and bottom stresses. To this end, an in-depth understanding of the near-bottom hydrodynamics and turbulence transition in the bottom boundary layer (BBL) under such waves is imperative, with the study of BBL instability being a prerequisite.

The solitary wave has been generally adopted as a canonical model of transient long waves in both experimental and theoretical studies since the 1970s. Although this idealization is not exactly representative of actual tsunami waves,^{3} Chan and Liu^{4} showed that the leading tsunami wave can be adequately modeled by a linear superposition of a number of *sech*^{2}-wave forms. In this study, as an additional simplification, we will restrict our analysis to the simple case where a single solitary wave is used to approximate small amplitude long waves typical of the coastal ocean.

Liu and Orfilla^{5} used a perturbation approach and the Boussinesq approximation to derive an analytical solution for the laminar viscous boundary layer flows under transient waves. Their analytical solution was based on the linearized boundary layer equations and can be found in the Appendix. Subsequently, Liu, Park, and Cowen^{6,7} extended this work to take into account the nonlinearity in the boundary layer equations and compared both analytical solutions, linear and nonlinear, with experimental data. They used particle image velocimetry (PIV) to measure the boundary layer velocity field under a solitary wave in a wave tank. They found good agreement between the theoretical solution and the experimental data. Another significant finding was that the boundary layer flow reverses its direction, and so does the bed shear stress, in the deceleration phase of the free stream velocity which always moves in the same direction of the wave propagation. The flow reversal is a result of the adverse pressure gradient which is established in the deceleration phase after the passing of wave crest. This pressure gradient gives rise to an inflection point in the velocity profile, which is a necessary but not sufficient condition for inviscid instability.^{8} In Liu *et al.*^{6} experiments, the free-stream Reynolds number, *Re*, of the experiments examined was *O*(10^{4}) and the boundary layer flow was always laminar. Specifically, the free-stream Reynolds number is defined as

where $ U o m \u2217 $ is the maximum wave-induced free-stream velocity, *ν*^{∗} is the fluid viscosity, and *a*^{∗} is the half of the wave-induced water particle displacement in the free-stream region, which is the ratio of maximum free-stream velocity $ U o m \u2217 $, to the wave angular frequency $ \omega \u2217 $ (to be defined in Sec. II).

Sumer *et al.*^{9} (hereafter SU10) pursued the experimental investigation of the boundary layer evolution under solitary waves with higher Reynolds number. To overcome laboratory limitations in attainable *Re*, they mimicked the BBL flow under a solitary wave using a U-shaped water tunnel in which the flow was driven by a time-varying, spatially independent, pressure gradient induced by a soliton-like wave in the flow direction. In this regard, it is emphasized that no actual surface solitary wave is present in the particular experimental study. Using the one-dimensional transient flow to simulate the inviscid free stream velocity under long waves is a reasonable approximation because, for very long waves, the spatial variation of the pressure gradient is small and can be neglected. Using this approach, SU10 were able to increase their Reynolds number to a maximum value of *Re* = 2 × 10^{6} which is two orders of magnitude higher than the previously discussed studies. SU10 showed that the boundary layer flow characteristics changed with increasing Reynolds number from laminar $ Re \u2009 < 2 \xd7 1 0 5 $ to laminar with two-dimensional vortex tubes near the bed $ 2 \xd7 1 0 5 \u2264 Re < 5 \xd7 1 0 5 $, and finally, to a transitional regime characterized by the appearance of single/multiple turbulent spots $ Re > 5 \xd7 1 0 5 $. However, due to the limitation of the experimental facility, SU10 were not able to reach a Reynolds number that is high enough for the BBL to become fully turbulent. SU10 demonstrated that their data agree well with the analytical solution of Liu and Orfilla^{5} and the previous wave tank measurements of Liu *et al.*^{6} for the laminar *Re* range (*Re* < 2 × 10^{5}). In the transitional regime, *Re* > 5 × 10^{5}, SU10 reported the appearance of turbulent spots, which appeared as single or multiple spikes in the bed shear stress measurements concurrently with the vortex tubes. SU10 speculated that at higher *Re*, i.e., *Re* > 2 × 10^{6}, the transition to turbulence can take place in the wave acceleration phase, although they could not explicitly reproduce it due to facility limitations. The mechanism behind the generation of the vortex tubes was interpreted by SU10 as the result of the inflectional-point shear layer instability similar to the observations of Carstensen *et al.*^{10} in oscillatory boundary-layer flows. It is worth noting that in their experiments, SU10 did not introduce explicitly any perturbation into their system but rather relied on the naturally available ambient disturbances in the laboratory environment.

Vittori and Blondeaux^{11} performed a three-dimensional direct numerical simulation (DNS) to simulate the boundary layer under a solitary wave. They discussed the dependence of the nature of the transition to turbulence on the ratio of wave amplitude to water depth. Following the SU10 study, Vittori and Blondeaux^{12} used the same numerical tool to investigate the development of coherent structures in the boundary layer and the efficacy of the log-law of the wall in approximating the structure of the near-bed flow field for a range of wave amplitudes. They reported a qualitative agreement between their numerical results and the experimental observations of SU10. They suggested that the generated vortices result from the growth of Tollmien-Schlichting waves (T-S). In each of the two studies,^{11,12} the base flow was perturbed by inserting weak random noise, of a relative order of amplitude *O*(10^{−4}) with respect to $ U o m \u2217 $, into the initially quiescent velocity field. Both studies were unable to simulate the formation of the turbulent spots that appeared in the BBL of the SU10 experiment. The lower bound for the emergence of the coherent structures was estimated to be somewhere below *Re* ≲ 5 × 10^{5}. Moreover, the transition to turbulence was suggested to take place immediately above this cutoff, i.e., *Re* > 5 × 10^{5}.^{12} They claimed that the discrepancy with the laboratory was due to the different sources of instability, i.e., bed roughness and vibrations, in the laboratory experiment in SU10. Furthermore, they indicated that the transition to turbulence and the formation of the coherent structures always took place in the wave deceleration phase.^{11,12} Scandura^{13} studied numerically the characteristics of the two-dimensional vortex tubes generated under solitary waves within a range of *Re* equivalent to SU10 experimental range. Scandura used small-amplitude sinusoidally shaped imperfections in the bottom boundary with a wavelength equivalent to the wavelength of the most unstable mode predicted from the linear stability analysis, as a source of disturbances. He showed that the mean spacing between the numerically developed vortex tubes is close to the one observed by SU10 for the range of *Re* and it was approximately constant for the range of the tested *Re*. The effect of the bottom wall perturbation is outside the scope of the investigation performed hereafter.

Ozdemir *et al.*^{14} carried out a number of three-dimensional direct numerical simulations to simulate the transition to turbulence in the BBL under a solitary wave for *Re* one order of magnitude higher than SU10. They inserted a two-dimensional random perturbation noise at the initial quiescent flow field with amplitudes varying between 1% and 20% of the maximum free stream velocity. Ozdemir *et al.* observed two possible instability mechanisms within a wave event: a short-lived non-linear long-wave instability that takes place during the acceleration phase for a large disturbance amplitude and a broadband linear instability which always takes place after the flow reversal during the wave decelerating phase. Moreover, they discussed the path of transition to nonlinear chaotic motions, i.e., turbulence, as a function of *Re*. They used a high order numerical model. Nevertheless, the transition to turbulence was always found to follow the breakdown of the span-wise rollers, i.e., vortex tubes, that results from the growth of the T-S in the deceleration phase. They divided the BBL flow under a solitary wave into four regimes: laminar, disturbed laminar, transitional and turbulent.

The observations of SU10 and previously discussed numerical studies indicate that an aspect of the wave-induced BBL which merits further investigation are its instability properties, namely, in the context of BBL response to externally imposed perturbations. The fundamental feature of these instability properties is the growth of 2D T-S waves which can be quantified, asymptotically, by solving the classical Orr-Sommerfeld (OS) equations for a single base flow profile. From the OS equations, the degree of perturbation amplification can be determined, either temporally or spatially, for steady parallel/nearly parallel base flow velocity profile. A classical example of the application of such a temporal stability analysis to a BBL flow is the case of the zero-pressure-gradient flat plate boundary layer.^{15} Following SU10’s assumptions, the BBL flow under solitary waves can be approximated, as previously indicated, as a sequence of temporally varying velocity profiles that are uniform in space, i.e., x-independent. This spatial uniformity satisfies the condition of the parallel base flow which allows the application of a temporal growth analysis focused on the OS equations instead of some other approach such as parabolic stability analysis.^{16} Nevertheless, as further elaborated in Sec. II, the unsteadiness of the BBL base flow must to be accounted for in the stability analysis. To this end, two approaches are commonly used: Floquet theory and the quasi-steady approximation. The former is only appropriate for periodic base flows as previously done for a large number of relevant studies.^{17–19} The latter, first introduced by Shen,^{20} relies on the assumption that the perturbation growth rate be faster than the rate of change of the base flow.^{21,22}

Blondeaux *et al.*^{23} were the first to carry out a quasi-steady linear stability analysis of the BBL flow at the bottom of a solitary wave using a “momentary” criterion of instability founded on the previously mentioned assumption of a perturbation development time scale that is much smaller than the time scale of the wave. Their analysis indicates that the laminar flow is linearly unstable for a range of flow parameters smaller than those observed in the SU10 (i.e., in the laminar regime, *Re* < 2 × 10^{5}). They suggested using the growth of the dimensionless kinetic energy per unit area of the flow perturbations as a measure for the flow transition from one regime to another (i.e., appearance of coherent structures). In their study, Blondeaux *et al.*^{23} reported fair agreement with the SU10 experiment in terms of the most unstable mode. However, the wave phase value at which the coherent structure appears in the experiment is much larger than the critical wave phase predicted by their linear stability analysis. We remark here that this difference could be caused by the fact that the perturbation amplitude needs some time to grow and reach a significant level in order to be observed, which contradicts the fundamental assumption of large separation between the time scale of the wave and that of the perturbation growth. In all these studies,^{11,12,23} instability is always observed in the deceleration phase for the range of wave parameters examined.

Recently, Verschaeve and Pedersen^{24} performed spatial stability analysis of the boundary layer under surface solitary waves by solving the parabolic stability equation. This technique allowed them to take into account the non-parallel effect of the base flow when including the non-linear term in the boundary layer equations. They also compare those results to the quasi-steady profile analysis, i.e., the Orr-Sommerfeld equation. Their preliminary results suggest that the approximation made by all previously referenced studies, i.e., neglecting the non-linear effect and studying the problem from the temporal point of view, is valid for the range of *Re* where the instability is significant. Investigating this problem in the spatial framework can be challenging in the sense of providing an adequate resolution of the different modes of instability.

The objective of this paper is to establish a more detailed map of instability in the flow parameter space for the BBL flow driven by a solitary wave. For the purpose of comparison with the laboratory observations of SU10, the base flow will be approximated according to the corresponding U-tube experiment: the transient base flow is driven by a uniform soliton-like wave induced pressure gradient.

The first step towards establishing such a map is to assess the limits of applicability of the quasi-steady linear stability analysis proposed by Blondeaux *et al.*,^{23} taking into account the time-dependent base flow. The next step consists of studying the variation in the resulting instability response as a function of the inserted noise characteristics (i.e., initial amplitude and insertion phase) for increasing *Re* using a fully non-linear two-dimensional pseudo-spectral code, which neither introduces nor amplifies any numerical noise. Additionally, the difference between the predicted instability in the theoretical/numerical studies and the one observed in the SU10 experiment is discussed in detail. Only two-dimensional analysis is performed. This allows one to reach a maximum *Re* higher than previous studies and to complete the stability map for the onset of two-dimensional instability.

Ultimately, in three dimensions, it is this two-dimensional instability which will lead to turbulence through the classical route of transition, at least according to Squires theorem.^{15} Note that there may exist different pathways towards transition, such as the formation of turbulent spots, which the two-dimensional approach followed here in this study will not be able to address.

The rest of the paper is structured as follows. Section II provides the problem formulation and governing equations. In Sec. III, a brief description of different numerical approaches used in the discretization of the governing equations is given. In Sec. IV, sample results are presented and then followed by a discussion in Sec. V. The conclusion is given in Sec. VI.

## II. PROBLEM FORMULATION

### A. Governing equations in the inviscid region

Consider a two-dimensional solitary wave of height *H*^{∗} that propagates in a constant water depth *h*^{∗}. Hereafter, the $ \u2217 $ sign indicates dimensional quantity. The Cartesian coordinate system $ x \u2217 , z \u2217 $ is employed, where *x*^{∗} points in the wave propagation direction and *z*^{∗} is the vertical coordinate pointing upward with its origin at the bed. Under such a wave, the irrotational velocity field, outside the boundary layer can be determined from potential flow theory.^{25}

As indicated earlier in the Introduction, this study will replicate the flow used in the U-tube experiment of SU10 as an approximation to the irrotational velocity component under a surface solitary wave. In this experiment, the generated free-stream velocity corresponds to the leading order solutions of Grimshaw’s formulae under the assumption of a small amplitude wave, i.e., $ H \u2217 / h \u2217 \u2192 0 $, and neglecting the relatively small spatial variation at any location under a long wave. Accordingly, the unidirectional horizontal velocity field outside the boundary layer can be expressed as follows:

where $ U o m \u2217 $ is the wave amplitude defined as

and *g*^{∗} is the gravitational acceleration. The wave frequency *ω*^{∗} is defined as

and, hence, a wave period (*T*^{∗}) can be defined as

which can be considered as the characteristic time scale of the solitary wave event.

Hereafter, we will omit the $ \u2217 $ from the product $ \omega \u2217 t \u2217 $ since this product is dimensionless. Strictly speaking, a solitary wave starts at *ωt* = − ∞ and ends at *ωt* = + ∞. In the following discussion, we shall truncate the entire event to −180^{∘} ≤ *ωt* ≤ 180^{∘}. The free stream starts from rest at *ωt* = − 180^{∘} and keeps accelerating until the maximum free stream velocity $ U o m \u2217 $ is reached, at the crest, i.e., *ωt* = 0. This acceleration stage is followed by a deceleration phase where the free stream velocity is reduced back to effectively zero due to an adverse pressure gradient.

Outside the BBL, the stream-wise pressure gradient from the x-momentum equation is simply

where *ρ*^{∗} is the fluid density. Fig. 1 shows a sketch of the free stream velocity normalized by the maximum free stream velocity $ U o m \u2217 $ and the corresponding pressure gradient normalized by the absolute maximum pressure gradient under a solitary wave.

### B. Governing equations in the boundary layer

#### 1. Fully nonlinear equations

The Navier-Stokes (NS) equations describing the motion of a three-dimensional, incompressible flow in the BBL in vector notation are

where $ u \u2192 \u2217 $ represents the velocity vector and *ν*^{∗} is the kinematic viscosity. In a two-dimensional domain (*x*^{∗}, *z*^{∗}), the velocity vector has two components: the horizontal stream-wise and vertical (*u*^{∗}, *w*^{∗}). The gradient vector is $\u2207= ( \u2202 \u2202 x \u2217 , \u2202 \u2202 z \u2217 ) $.

Consequently, the total pressure (*p*^{∗}) in Eq. (7) can be divided into two parts,

where *p*^{″∗} is the dynamic pressure that develops in response to the BBL flow and later to the velocity perturbation. As in the U-tube experiment, carried out by SU10, the flow is driven by a pressure gradient $\u2202 p 0 \u2217 /\u2202 x \u2217 $ given in Eq. (6) that is constant in space yet variable in time.

To nondimensionalize the problem, we use $ U o m \u2217 $ as a velocity scale and the Stokes laminar boundary layer thickness *δ*^{∗},

as a length scale.

Using the above length and velocity scales, the non-dimensional NS equations in the BBL in Cartesian coordinates are written as

where the non-dimensional quantities are non-starred and the BBL Reynolds number based on the boundary layer thickness is

We can express the relation between the BBL Reynolds number *Re*_{δ} and the free stream velocity based Reynolds number *Re*, defined in Eq. (1), as

#### 2. Linear stability analysis formulation

Linear stability analysis can be employed to study the stability of a BBL flow under the assumption that the base flow is parallel or nearly parallel. A traditional quasi-steady linear stability analysis, similar to the one by Blondeaux *et al.*,^{23} is first considered in this paper where a “momentary” criterion for instability is assumed and an eigenvalue problem is solved. Hence, in order to take into account the transient nature of the base flow, especially for the relatively low *Re*_{δ} regime, an initial value problem is solved in place of the eigenvalue problem.

##### a. Quasi-steady analysis: Eigenvalue problem.

For the sake of completeness, a brief review of the derivation of the OS equation (see Schmid and Henningson^{26} for further details) is presented hereafter for a steady velocity profile. The flow field is decomposed into a laminar steady base flow and a fluctuating component as follows:

In the above Eqs. (16) and (17), the time dependence of the base flow is neglected and *P*(*x*) and *U*(*z*) represent the quasi-steady base flow mean pressure and velocity. The unperturbed horizontal velocity in the boundary layer for a specific *ωt*, i.e., *U*(*z*), is calculated from the analytical solution of Liu and Orfilla.^{5} The tilde indicates the fluctuating component and the uppercase letter represents the base flow. Substituting Eqs. (16) and (17) in Eqs. (11) and (12), subtracting the leading order balance of the base flow and neglecting the non-linear terms, results in the following linearized equations:

The (′) denotes derivative with respect to *z*. Eliminating the pressure terms in the above equations, the equation for the stream-wise vertical velocity ($ w \u0303 $) is obtained as follows:

The boundary conditions $ w \u0303 =\u2202 w \u0303 /\u2202z=0$ both at the bed $ z = 0 $ and in the far-field $ z \u2192 \u221e $ are required. In this study, the far field upper boundary of the domain of Eq. (21) is fixed at *z* = 15 (dimensionless), which is sufficiently far from the bed and at this location, the flow can be reasonably assumed to be uniform in the vertical direction. To carry out the quasi-steady analysis according to the classical modal approach,^{26} the perturbation velocity is decomposed into normal modes as follows:

where *α* is the wavenumber in the *x* direction ($\alpha = 2 \pi \lambda x ,\u2009 \lambda x $ is the stream-wise wavelength of the perturbation, which eventually determines the spacing between the resulting vortex tubes), $ w \xaf $ is the perturbation wave-like amplitude, and *ω _{p}* is the perturbation frequency $ \omega p = \omega r + i \omega i $. Substituting Eq. (22) in Eq. (21) (equivalent to taking the Fourier transform in the horizontal directions), the well-known OS equation is obtained as follows:

with the OS operator *L _{OS}* defined as

In the equation above, *D _{z}* is the derivative operator, with respect to

*z*.

*k*=

*α*for a two-dimensional case. For the temporal stability analysis,

*ω*is complex and it appears as the eigenvalue in the OS equation for each profile. If

_{p}*ω*> 0, this indicates an unstable mode (i.e., growth in time), whereas

_{i}*ω*< 0 corresponds to a decaying mode.

_{i}Five important quantities are used to characterize the instability properties of a particular *Re*_{δ} as estimated from quasi-steady linear stability analysis. The earliest time at which growth contours of *ω _{i}* cross the zero threshold towards monotonically increasing positive values is regarded as the critical wave phase, $\omega t critical $ which has its corresponding critical wavenumber $ \alpha critical $. The maximum growth rate,

*ω*

_{imax}, occurs at $\omega t \omega i max $ and for a corresponding wavenumber $\alpha \omega i max $.

Recall from the Introduction that the main assumption of using the quasi-steady analysis for a transient base flow is that the perturbation growth rate is faster than the rate of change of the base flow. As will be shown in Sec. IV A, the applicability of this assumption will be tested for the particular BBL of interest by comparing both quasi-steady and transient analyses.

##### b. Transient analysis: Initial value problem.

If the transient behaviour of the base flow must be taken into account, an initial value problem must be solved where a different solution form of the perturbation velocity is sought for

the linear initial value problem equation (26), which can be re-written as

where $ L 1 =\u2212 L O S / k 2 \u2212 D z 2 $. It must be emphasized that this transient analysis can be only performed for one mode at a time. To monitor the energy growth of the perturbation with time, the vertical perturbation kinetic energy *q* is introduced and, for this case, can be calculated from the integration of $ w \u0302 ( z , t ) $ over the domain volume *V* as follows:

Quasi-steady analysis indicates that positive *ω _{i}* never occurs before the flow reversal phase at the range of lower

*Re*

_{δ}where transient stability analysis is necessary (Sec. IV A). Therefore, when performing this particular analysis,we choose to perturb the base flow starting from the wave crest phase, i.e.,

*ωt*= 0

^{∘}, with a maximum relative initial amplitude of $O 1 0 \u2212 3 $.

## III. NUMERICAL METHOD

Since the presentation of our results first explores the applicability of the quasi-steady linear stability analysis as a function of Reynolds number, the corresponding numerical tools are considered first here, followed by an overview of the computational scheme used within the fully nonlinear flow solver.

### A. Linear stability analysis

For the numerical solution of either the classical eigenvalue problem, Eq. (23), or the initial value problem, Eq. (26), a Chebyshev collocation method^{27} is used in the discretization of the linear operators of the vertical velocity perturbation in Eq. (21), i.e., $ D z 2 $ and $ D z 4 $, on Gauss-Lobatto-Chebyshev grid points (GLC) as recommended by Orszag^{28} in order to ensure high accuracy.

#### 1. Quasi-steady analysis: Eigenvalue problem

After assuming a wavelike solution for the perturbation, according to Eq. (22), Orr-Sommerfeld equation (23) is solved to find the asymptotic solution for the individual frozen base velocity profiles under the “momentary” condition of instability proposed by Blondeaux *et al.*^{23} The solution of such an equation consists of computing *ω _{i}* > 0 according to (23),

This equation is solved for each individual base flow velocity profile, assuming that it is steady, and supplies a growth rate, *ω _{i}*, as a function of the perturbation wavenumber

*α*.

#### 2. Transient analysis: Initial value problem

Alternatively, for the initial value problem, the evolution of the base flow is first sampled over equally spaced intervals of small duration $\Delta \omega t =0. 1 \u2218 $ within which it is safe to assume that the base flow is nearly steady. The length of this interval, i.e., $\Delta \omega t $, can be translated into an equivalent perturbation growth time, i.e., *t _{p}*, using the problem time scale $ i.e., \u2009 \delta \u2217 / U 0 m $. Within this $\Delta \omega t $ interval, the linear operator

*L*

_{1}is considered constant for an individual mode. The Crank-Nicolson scheme is then used in the temporal discretization of Eq. (27), using an equal time step of Δ

*t*to proceed from step

_{p}*n*to step

*n*+ 1 for the each of the total interval time

*t*as follows:

_{p}This approach takes into account the effect of the transient nature of the base flow on the cumulative growth of the perturbation. Recall that this analysis is only performed for a single streamwise mode. However, in this type of analysis, the transient growth is highly dependent on the vertical structure of the initially inserted perturbation fields and on how close it may be to the most unstable eigenvector. To eliminate any resulting ambiguity, the evolution of the perturbation that is ultimately reported is the ensemble average over a total number of *N* simulations each of which has a different initial perturbation field with a random spatial distribution.

### B. Fully non-linear simulation

Two-dimensional fully non-linear governing equations (11) and (12) are solved in a computational domain equivalent to an idealized box, corresponding to a fixed observer under the wave, where the velocity field is driven by a forcing term given by Eq. (6). The lateral domain boundaries are periodic boundary conditions, whereas the top and the bottom boundaries are free-slip and no-slip boundary conditions, respectively. The resulting velocity field is then perturbed at an appropriately chosen time of insertion with a prescribed numerical noise field.

The numerical method used to solve the problem governing equations is a spectral multidomain penalty method model.^{29} This pseudo-spectral model has been used to examine similar boundary layer flows associated with internal solitary waves,^{30} albeit in a rather different problem geometry. More details regarding the numerical model may be found in Diamessis *et al.*^{29} The semi-implicit temporal discretization relies on a third-order stiffly stable scheme.^{31} Additionally, an adaptive time stepping scheme is employed to decrease the time step size during the more energetic flow phases when the instability kicks in.

The spatial discretization relies on a Fourier discretization in the horizontal direction, where *N _{x}* is the total number of equally spaced grid nodes in the spanwise direction, and on a Legendre-polynomial-based single/multidomain discretization in the vertical. In the majority of the cases simulated throughout this study, only a single vertical domain, i.e., no subdomains, is employed to avoid any numerical instability that may be introduced due to discontinuity at the subdomain interface in the highest Re cases considered. Additionally, explicit spectral filtering is employed in both spatial directions to eliminate any spurious numerical noise and preserve the high accuracy of the model.

^{32}In all simulations, the horizontal grid spacing, Δ

*x*, is less than

*δ*/3. In the vertical direction, at least between 12 and 15 grid nodes are located within one Stokes thickness

*δ*to enable sufficient resolution of the developing BBL and of the subsequently generated vortex tubes.

For the purpose of validation, the numerical results of laminar cases, i.e., simulations in the absence of any external perturbation, are in a very good agreement with the analytical solution (see the Appendix).

## IV. RESULTS

The presentation of our results is structured by our intention to illustrate the instability properties and associated dynamical response of the transient BBL under examination as a function of increasing *Re*_{δ}, much in alignment with the approach of SU10. First, the limits of applicability of the quasi-steady linear stability analysis, as a mean of interpretation of the perturbation growth characterisers in the linear regime, are explored in the context of the transient nature of the base flow. The fully nonlinear solver is then used to examine the non-linear phase of BBL evolution which consists of the formation of vortex tubes and any interactions among them. Having demonstrated agreement with SU10’s findings, we extend our study to a higher *Re*_{δ}, with a focus on the connection between the spacing of the observed vortex tubes and the most unstable mode as established through both linear stability analysis and fully nonlinear simulations.

### A. Applicability of the quasi-steady linear stability analysis

In their study, Blondeaux *et al.*^{23} suggest that for small wave non-linearity (corresponding to small values of *H*^{∗}/*h*^{∗}, which they indicate to be the critical parameter determining the stability characteristics of the solitary wave-induced BBL), the “momentary” criterion of instability is not applicable. Hence, in this regime, quasi-steady linear stability analysis is not expected to provide accurate results. However, the precise upper bound of this regime, expressed in terms of *H*^{∗}/*h*^{∗}, was not specified by Blondeaux *et al.*^{23}

In this regard, the scaling introduced in Eqs. (12)-(15) indicates that, at least for the transient BBL examined in this paper, *Re*_{δ} is the sole non-dimensional parameter governing its instability response. To illustrate the presence of an upper bound of *Re*_{δ}, below which quasi-steady linear stability analysis is not applicable, the case with *Re*_{δ} = 211 (equivalent to *Re* = 2 × 10^{4}) is now examined with an emphasis on the disagreement between the quasi-steady and time-dependent approaches. The applicability criterion is to be formally defined in the next page.

Quasi-steady asymptotic stability analysis is used to construct the contour map for *ω _{i}* > 0 throughout the wave deceleration phase as shown in Fig. 2. Recalling the definitions in Sec. II B 2, the maximum growth rate is

*ω*

_{imax}= 0.0128 and the critical wave phase occurs at

*ωt*≈ 40

_{c}^{∘}. The critical wave phase happens after the flow reversal which takes place at

*ωt*≈ 30

^{∘}for all cases.

The next step is to solve the initial value problem associated with the transient analysis for the same case, i.e, *Re*_{δ} = 211. As pointed out before, this type of analysis is only performed for an individual mode at once. Hence, we shall show the results for the first destabilized mode, i.e., the critical wavenumber, $ \alpha critical =0.44$, as identified by the quasi-steady analysis. Across all *Re*_{δ}, the highest growth rate modes, i.e., $\alpha \omega i max $, always lies within the range of $\alpha \u2208 0 . 43 , 0 . 47 $. The critical wavenumber for *Re*_{δ} = 211 happens to also be in this range as shown in Fig. 2.

The default approach for the transient stability analysis is to insert the perturbation at *ωt* = 0^{∘} since any growth is only anticipated during the wave deceleration phase for this relatively low *Re*_{δ} value. The amplification factor, *G*, is examined as a measure of the growth of the inserted perturbation amplitude. *G* is defined as the kinetic energy based on the vertical velocity component, $q t $, normalized by the initial energy level, *q*_{0}, at the time of insertion.

The evolution of the amplification factor for three modes representative of the highest growth rate *α*-interval defined above is shown in Fig. 3. The inserted perturbation amplitude undergoes an initial steep decay $ 0 \u2218 \u2264 \omega t \u2264 4 \u2218 $, and is then followed by a mild decay up to *ωt* = 100^{∘} and then remains nearly constant. Although the quasi-steady analysis indicates a considerable perturbation growth, the transient analysis shows no growth but instead a continuous decay. Results from the corresponding two-dimensional DNS for this case, initialized with a low-amplitude white noise spectrum (not shown here), confirm this continuous decay.

Having illustrated the disagreement between quasi-steady and transient stability analyses at a particular value of sufficiently low *Re*_{δ}, the lower limit of applicability of the former as a reliable descriptor of the instability of the transient BBL flow is identified through a series of numerical experiments over a range of gradually increasing *Re*_{δ}. The following applicability criterion is specifically proposed: the quasi-steady analysis is deemed to be valid if the degree of perturbation growth is such that, at least by the end of the wave event, the perturbation amplitude exceeds its original value. In other words, the amplification factor of the inserted perturbation has to become equal to one within the period of the solitary wave. Nevertheless, satisfaction of such a criterion does not guarantee the formation/appearance of vortex tubes as will be elaborated in Sec. IV B.

A representative example of amplification factor *G* evolution curves for different *Re*_{δ} is shown in Fig. 4. The proposed criterion is first satisfied at *Re*_{δ} = 350. The growth rate contour for this *Re*_{δ} case is shown in Fig. 5. What is interesting in this last figure is that the critical wave phase is $\omega t c r i t i c a l \u22483 0 \u2218 $ when the flow reversal takes place. Therefore, a consequence of the amplification factor based criterion is that the applicability of the quasi-steady linear stability analysis for the transient base flow is intimately linked to the neutral curve crosses the near-bed flow reversal development threshold, i.e., *ωt* ≈ 30^{∘}. The applicability limit defined earlier is confirmed using DNS, which shows trends equivalent to those in Fig. 4. In DNS, the amplification factor of different individual modes are monitored. The modes are isolated by Fast Fourier transformation (FFT), of the data in the streamwise direction at every time.

The next question is whether the T-S waves of instability can grow adequately to be noticeable, i.e., vortex tubes actually form. This happens when the amplitude of the velocity perturbation reaches the same order of magnitude of the base flow. In this case, the non-linearity will take over and will re-distribute non-negligible energy across different modes giving rise to the formation of vortex tubes. Vortex tube formation is, of course, hard to identify using any of the linear stability analysis alone since we neglect the non-linear term [Eqs. (19) and (20)] and, consequently, any kind of interaction among individual modes. The correct reproduction of nonlinear effects motivates the use of DNS to further explore the parameter range identified as linearly unstable.

### B. Two-dimensional direct numerical simulation analysis

In this analysis, the goal is to identify the *Re*_{δ} limit above which the velocity perturbation will grow to an order of magnitude comparable to the base flow. The base flow then deviates from its laminar state as a result of the formation of coherent vortex structures. To trigger the instability in the numerical simulations, white noise is inserted with a prescribed small amplitude. This approach allows the reproduction of a wide range of possible instability modes. The inserted noise is designed to satisfy the divergence-free condition which already holds for the base flow. A sensitivity analysis (not shown here) has indicated that the two most important parameters associated with the introduced noise are its initial amplitude and insertion phase.

It is important to note that the initial amplitude of the inserted perturbation has to be much smaller than that of the base flow at the time of insertion. However, this was not the case in some previous studies.^{12,14} Use of a high initial amplitude noise, as compared to the base flow at the time of insertion, can induce some premature non-linear effects and can modify the linear instability response of the base flow in the wave deceleration phase. In this regard, the experiments of SU10, where presumably any ambient perturbations are relatively weak, did not report any deviation from the laminar state throughout the acceleration phase for their *Re*_{δ} value range, as indicated by their time series of the bed shear stress.

On the other hand, in the DNS, the insertion phase of the perturbation can significantly impact the time of appearance of the vortex tubes. For instance, when a small perturbation is inserted at the wave initiation phase, i.e., *ωt* = − 180^{∘}, the perturbation amplitude will undergo a steady viscous decay throughout the acceleration phase, namely, for relatively lower *Re*_{δ} values. In this case, the perturbation amplitude at the critical wave phase is very small. Although the perturbation will eventually grow after this point, it will not be able to reach the order of magnitude of the base flow within the time scale of the wave event and the vortex tubes will never form.

Similarly, SU10 adopted the deviation of base flow from its laminar state and the associated flow structure to serve as a basis for the classification of the BBL flow under solitary wave. Such a deviation results from the base flow interaction with the same order of magnitude perturbation. This particular, more qualitative, classification is motivated by the capacity to observe the manifestation of the well-developed perturbation, in the form of generated vortex tubes, instead of identifying instability more quantitatively, as done using linear stability analysis through its focus on perturbation growth rates (Sec. IV A). Accordingly, in the discussion that follows, the linearly unstable regime of flow behaviour, earlier identified from the transient stability analysis, will be redivided into three regimes:

Stable (S) or Laminar: no vortex tubes are formed regardless of the perturbation growth.

Conditionally Unstable (CU): the formation of the vortex tubes within the wave event is conditional on the initial perturbation amplitude and insertion phase.

Unconditionally Unstable (UU): vortex tubes are formed regardless of the perturbation characteristics (i.e., initial amplitude or insertion time).

#### 1. Conditionally unstable regime (450 < *Re*_{δ} < 1900)

During the initial stage of this study, a general destabilization technique similar to what was adopted in some of the previous numerical studies (e.g., Vittori and Blondeaux^{12} and Ozdemir *et al.*^{14}) was considered.

As elaborated in Sec. IV B, no perturbations are inserted at the very initial wave phase, i.e., *ωt* = − 180^{∘}. When low in amplitude, these perturbations simply decay due to viscosity and cannot trigger any instability in the deceleration phase, whereas, when high, these perturbations can induce unrealistic non-linear instability during the acceleration phase that may have some effect on the subsequent deceleration phase instability as reported in a previous study.^{14}

Instead, the perturbation with the desired amplitude is inserted at the wave crest, i.e., *ωt* = 0^{∘}, as is done in the transient linear stability analysis (Sec. IV A). The extensive sensitivity analysis performed (not shown) has also indicated that inserting the noise at the initial stage or at the crest will result in the same response provided that the base flow is stable from the linear stability point of view. Additional numerical experimentation on our side has shown that the appropriate maximum relative initial amplitude of the inserted perturbation is $O 1 0 \u2212 3 $ to avoid the non-linear modification of the base flow.

Experimental results and observations at *Re*_{δ} = 920 (equivalent to test number 8 in SU10, *Re* = 4.4 × 10^{5}) provide a benchmark for comparison with numerical results for the conditionally unstable regime. For this specific *Re*_{δ} value, SU10 report the formation of an array of vortex tubes late in the deceleration phase. They also provide the time series of the bed shear stress measurement of this test case. In the corresponding numerical simulation, a maximum initial perturbation amplitude of $O 1 0 \u2212 3 $ was inserted at the wave crest phase. The formation of vortex tubes near the bed was observed at *ωt* = 80^{∘}. The distance between two successive vortex tubes, which originate from the most unstable mode, for this *Re*_{δ} value case is λ_{x} ≈ 14.5, which corresponds to *α _{x}* = 2

*π*/λ

*x*≈ 0.43. SU10 reported an average spacing between the vortex tubes λ

_{sumer}≈ 14 for this same range of

*Re*

_{δ}which is in good agreement with the distance obtained from our numerical simulations.

From the linear stability point of view, *α _{x}* = 2

*π*/λ

*x*≈ 0.43 is in close proximity to both the critical wavenumber and the one having the largest growth rate (Fig. 6; see also the discussion on highest growth rate modes interval in Sec. IV A). An even more rigorous comparison is offered in Fig. 7, which shows the time series of both measured and calculated bed shear stresses. The DNS time series of bed shear stress, sampled at different locations along the bed, agrees well with the corresponding time series from the experiment. Moreover, the phase at which the bed shear stress starts to deviate from the laminar solution, i.e.,

*ωt*≈ 80

^{∘}, is related to the roll-up phase of the vortex tubes. For the evolution of both total vorticity and perturbation vorticity fields from the DNS simulation of

*Re*

_{δ}= 1100, as a representative sample for the conditionally unstable regime, see supplementary movies 1 and 2.

^{33}

The other point to address here is the determination of the lower limit of this instability regime where vortex tubes are observed in the deceleration phase. Moreover, the source of any discrepancy with the experiment is also pursued. SU10 identified the lower limit of this regime, i.e., laminar flow with vortex tubes, as $Re=2\xd71 0 5 \u2009 Re \delta = 630 $. In contrast, Vittori and Blondeaux^{11,12} identified a slightly higher Reynolds number limit, $Re=5\xd71 0 5 Re \delta = 1 0 3 $. Vittori and Blondeaux suggested that the reason for their higher limit was due to presence of additional external sources of perturbations in the experiment (i.e., vibration). Although one may partially agree with this explanation, an additional interpretation behind their lack of observation of vortex tubes at lower Re is due to the insertion of their initial *O*(10^{−4}) perturbation into the initially quiescent velocity field, at *ωt* = − 180^{∘}. Hence, the perturbation undergoes a strong viscous decay in the acceleration phase as explained earlier. This interpretation has been confirmed by the subsequent work of Ozdemir *et al.*^{14} who also inserted the perturbation at the initially quiescent velocity field and observed the vortex tubes only when the amplitude of the initial perturbation was relatively high (approximately 20% of the max free stream velocity *U*_{0m}).

As for the different *Re*_{δ} limit identified by SU10 compared to Vittori and Blondeaux,^{11,12} and apparently this study, a more rigorous explanation may be offered. Following the same approach adopted for *Re*_{δ} = 920, i.e., inserting the maximum possible perturbation amplitude $O 1 0 \u2212 3 $ at the crest, the time series of the DNS bed shear stress is plotted against the laminar solution (Liu *et al.*^{6}) for *Re*_{δ} = 450 case in Fig. 8. This figure is very similar to figure number 6 in SU10 where they compared the time series of the measured bed shear stress to the analytical laminar solution for $Re=9.4\xd71 0 4 \u2009 Re \delta = 430 $ case. SU10 disregarded their measurements/visualizations for *ωt* > 100^{∘} due to the experimental limitation particular to the generation mechanism of the base flow driven by a soliton-like wave induced pressure gradient. In Fig. 8, the bed shear stress indeed follows exactly the laminar solution up to *ωt* = 100^{∘}. This explains why this *Re* case was classified by SU10 to be laminar.

However, the distortion observed in the DNS bed shear stress for *ωt* > 100^{∘} reveals an interaction between the amplified perturbation and the base flow. This interaction manifests itself in a significant distortion of the near-bed parallel shear layer giving rise to the very initial stages of formation of a coherent vortex tube. According to our classification of instability regimes in the beginning of Sec. IV B, this *Re*_{δ} value is considered unstable. For any *Re*_{δ} ≳ 450, with $O 1 0 \u2212 3 $ perturbations introduced at *ωt* = 0^{∘}, vortex tube arrays are expected to form as the result of the breakdown of the shear layer within the time scale of the wave event.

#### 2. Unconditionally unstable regime (*Re*_{δ} > 1900)

In this regime, the perturbation is inserted into the initially quiescent velocity field of the wave, i.e., *ωt* = − 180^{∘}; viscous decay is not sufficiently strong to damp perturbation growth. Very low amplitude perturbation must be used in order to avoid the excitation of any non-linear instability in the acceleration phase. The non-dimensional relative initial perturbation amplitude is $O 1 0 \u2212 12 $ to $O 1 0 \u2212 8 $. To determine the lower *Re*_{δ} bound of this regime, a series of numerical runs were performed while increasing *Re*_{δ}. Since only the two-dimensional onset of instability and the formation of the coherent structures leading to the classical transition to turbulence are of interest, the base flow and perturbation are uniform in the span-wise direction. Hence, the formation of the turbulent spots, which can be an alternative mechanism to transition as observed in SU10, cannot be simulated.

Starting with *Re*_{δ} = 1900, the growth rate of the perturbation is high enough that the vortices will eventually evolve regardless of the initial amplitude of the perturbation as explained before. The unconditionally stable regime can be further divided into two main sub-regimes where the difference between them will be discussed later. The first sub-regime is for the interval 1900 < *Re*_{δ} < 5000, whereas the second sub-regime is for the interval *Re*_{δ} > 5000. In the first sub-regime, the dominant unstable perturbation mode, i.e., the mode that sets the spacing between the generated vortex tubes, was found to be also the mode with the highest growth rate as computed by the quasi-steady linear stability analysis.

This mode has a wavelength of λ_{x} ≈ 15 and it predominantly grows after the flow reversal. In the second sub-regime, i.e., *Re*_{δ} > 5000, the wavelength of the most unstable mode increases to λ_{x} ≈ 21. This wavelength corresponds to the first destabilized mode, although it is not necessarily the fastest growing one. This longer mode predominantly grows before the flow reversal phase (i.e., *ωt* < 30^{∘}). When *Re*_{δ} ≥ 7900, the base flow starts to be destabilized in the acceleration phase of the wave, before the occurrence of the adverse pressure gradient and the associated formation of the inflectional point in the velocity profile. Such an instability response suggests a viscous origin of the instability, similar to what is observed in the plane Poiseuille or Blasius boundary layer flows (see Criminale *et al.*^{8} for details).

The evolution of the amplification factor *G _{W}* of each of the most unstable modes for the different

*Re*

_{δ}cases is shown in Fig. 9.

*G*is the equivalent of the amplification factor

_{W}*G*used in the transient linear stability analysis (Sec. IV A), although it is defined for a particular mode of interest and not for the entire flow field. The dotted curves in Fig. 9 represent the transient variation of

*G*for dominant perturbation wavelength, i.e., λ

_{W}_{x}≈ 15, in the first sub-regime, i.e., 1900 <

*Re*

_{δ}< 5000. The solid curves are representative of the transient variation of

*G*in the second sub-regime, i.e.,

_{W}*Re*

_{δ}> 5000, where the wavelength of the most unstable mode is λ

_{x}≈ 21. As shown in this figure,

*G*undergoes a consistent growth until it reaches a plateau at the time when the coherent vortex tubes forms. The figure also shows that the growth rate of the most unstable mode increases as the base flow

_{W}*Re*

_{δ}increases, which agrees with the quasi-steady linear stability analysis (e.g., Figs. 5 and 6).

## V. DISCUSSION

### A. Mode of instability

As previously discussed, depending on *Re*_{δ}, there are two distinct dominant instability modes. The first mode, hereafter regarded as *M*1, is associated with the highest energy growth rate, which consistently takes place after the flow reversal. The wavenumber of this mode is roughly around $\alpha \u2208 0 . 47 , 0 . 42 $ which corresponds to a non-dimensional wavelength of $ \lambda x \u2208 13 . 5 , 15 $. Above a certain *Re*_{δ} threshold, another mode, hereafter regarded *M*2, is the first to destabilize before the flow reversal and has a wavenumber of roughly *α* ≈ 0.30 that corresponds to a non-dimensional wavelength of λ_{x} ≈ 21. The maximum growth rate of *M*2 is much smaller than the highest energy growth rate associated with *M*1 throughout the BBL evolution. However, as described in Sec. IV B 2 for *Re*_{δ} > 5000, M2 can grow sufficiently fast such that it prevails over *M*1 as the dominant instability mode.

To further highlight the difference between these two instability modes, the base flow neutral contour curves, i.e., *ω _{i}* = 0, are shown in Fig. 10 as a function of

*Re*

_{δ}following the calculations of quasi-steady linear stability analysis. In this figure, the locus of centers of the highest growth contour for the different

*Re*

_{δ}is represented as a solid ellipse. For the conditionally unstable regime 450 ≲

*Re*

_{δ}< 1900, the first destabilized mode is the one with the highest growth rate and corresponds to

*M*1 observed in the two-dimensional DNS (e.g., Fig. 9). This numerical observation of

*M*1 dominance fits well with the average distance between the successive vortex tubes observed in the SU10 study, i.e., $ \lambda Sumer \u2208 13 , 14 . 5 $, and in previous numerical simulations

^{12–14}within this particular

*Re*

_{δ}range. In the first sub-regime of the unconditionally unstable regime, i.e., 1900 <

*Re*

_{δ}< 5000, although

*M*2 is the first destabilized mode, it does not have enough time to grow to a significant amplitude before the flow reversal and therefore control the formation of the coherent vortex structures. In this

*Re*

_{δ}range, the

*M*1 mode, which always has the highest growth rate, is the one that sets the spacing between the vortex tubes formed after the flow reversal. This instability is equivalent to the one reported for the conditionally unstable regime. Consistently with our findings, Ozdemir

*et al.*

^{14}reported a comparable average distance between the formed vortex tubes, i.e. λ

_{Ozdemir}≈ 13.5 for

*Re*

_{δ}= 2000 and λ

_{Ozdemir}≈ 16 for

*Re*

_{δ}= 2500.

As *Re*_{δ} is increased further, the neutral curves start to show a visible protrusion towards a smaller wavenumber that of *M*2. Eventually, this protrusion penetrates into the acceleration phase region. As discussed earlier, for *Re*_{δ} > 5000, *M*2 starts to grow at an earlier phase and it can actually overtake the high growth rate of *M*1 that occurs after the flow reversal.

To further illustrate the competition between the two modes, two values of *Re*_{δ}, 4000 and 6300, are additionally examined. In Fig. 11, the growth contours for *Re*_{δ} = 4000 are first plotted. Fig. 12 shows the corresponding transient growth rate of *G _{W}* for the two modes

*M*1 and

*M*2. For this latter

*Re*

_{δ}value case,

*M*2 has the highest growth rate prior to the flow reversal, i.e.,

*ωt*< 30. After the flow reversal takes place,

*M*1 dominates, as induced by the increasingly steeper slope, and eventually, the generated vortex tubes are spaced in agreement with the quasi-steady instability theory-computed wavelength of

*M*1 (i.e., λ

_{x}≈ 15; see also the left panel of Fig. 15 and its discussion).

An additional exercise which further demonstrates the competition between the two instability modes is conducted at *Re*_{δ} = 6300 where the perturbation is intentionally inserted at three different wave phases. The growth rate contours of this *Re*_{δ} case, obtained from the linear stability analysis, are plotted in Fig. 13. Three vertical lines, representing the different insertion phases of perturbation: *ωt _{ins}* = 0

^{∘}, 18

^{∘}, and 30

^{∘}are also included in this figure.

At the same value of *Re*_{δ}, Fig. 14 shows the transient variation of *G _{W}* for each

*M*1 and

*M*2 for the three different times of perturbation insertion $ i.e. , \u2009 \omega t ins = 0 \u2218 , \u2009 1 8 \u2218 , \u2009 and \u2009 3 0 \u2218 $ where the initial relative amplitude of the inserted perturbation is $\u223cO 1 0 \u2212 12 $. In the calculation of

*G*for the different cases of insertion, the kinetic energy of the vertical velocity component for the particular mode (i.e., $ q W \omega t $) is normalized by its initial value at the insertion time $ q W \omega t ins $. In the left panel of Fig. 14,

_{W}*M*2 starts to grow immediately at

*ωt*= 0

^{∘}, as the adverse pressure gradient is established, with a relatively high growth rate compared to

*M*1.

*M*1 picks up around the flow reversal phase, i.e.,

*ωt*≈ 30

^{∘}, but cannot overcome the earlier growth of

*M*2 (as in the case

*Re*

_{δ}= 4000, Fig. 12) provided that the perturbation was inserted prior to flow reversal $ i.e., \u2009 \omega t \u2248 3 0 \u2218 $. As the perturbation is inserted closer to the flow reversal phase angle, the most unstable mode changes from

*M*2 to

*M*1 as shown for the

*ωt*≈ 30

_{ins}^{∘}case.

In Fig. 15, a comparison between two different *Re*_{δ} values, i.e., 4000 (Figs. 15(a)–15(c)) and 8000 (Figs. 15(d)–15(f)), in terms of the spatial structure of instability modes throughout the different development phases of the perturbation, is shown. Stream-depth slices of the perturbation vorticity contour, computed from both linear stability analysis (Figs. 15(a) and 15(d)) and DNS (Figs. 15(b), 15(c), 15(e), and 15(f)), are used to highlight the difference between those two *Re*_{δ} cases. In each of these contour plots, the vorticity field is normalized by the local maximum at the corresponding phase. Figs. 15(a) and 15(d) represent the vorticity field of the most unstable mode, at the indicated phases, estimated from the linear stability analysis.

The first two rows of Fig. 15 indicate a very good agreement between the linear stability analysis and the DNS results in terms of the structure of the perturbation vorticity field. From these first two rows, two shear layers exist on top of each other. The upper layer is similar to the classical *tanh*-profile shear layer, whereas the lower one is induced by the no-slip constraint at the wall. As *Re*_{δ} increases, the lower layer shape is much sharper, i.e., sheared, and compacted vertically as the perturbation starts to evolve at earlier times. At a later stage of the perturbation evolution, the spacing between the formed vortex tubes is set by the dominant instability as mentioned earlier. *M*1 dominates for *Re*_{δ} = 4000, whereas *M*2 is the dominant mode for *Re*_{δ} = 8000. The difference in spacing between the resulting vortex tubes can be observed in the two plots of the last row, i.e., Figs. 15(c) and 15(f) where the diameter of the vortex tubes is directly proportional to the thickness of the base flow shear layer at the rolling up phase of the vortices. Hence, for lower *Re*_{δ} and late emergence phase, the diameter of the vortex tubes is larger than at higher *Re*_{δ} case.

### B. Similarity with instability properties of other wave-driven transient flows

The *Re*_{δ} limit below which quasi-steady instability analysis is not applicable may be generalized to different wave-induced flows. In particular, this generalization may be made for flows driven by internal waves, albeit controlled by a different governing parameter. In this regard, both Troy and Koseff^{34} and Barad and Fringer^{35} reported an instability criterion based on the non-dimensional average perturbation growth rate for progressive internal gravity waves and for internal solitary waves, respectively. In both cases, the particular instability criterion is satisfied when $ \sigma i \u0304 T w >5$, where $ \sigma i \u0304 $ is the growth rate of the instabilities averaged over a time interval *T _{w}* during which a parcel of fluid within stratified water is subjected to a wave-driven background vertical shear corresponding to Richardson number

^{36}less than 1/4 and the base flow was assumed to be parallel as in the current study. For the problem at hand, a similar non-dimensional average growth rate $ \omega i \u0304 T W I $ may be calculated where $ \omega i \u0304 $ is average positive growth of the critical wavenumber, i.e., $ \alpha critical $, found to always lie within the highest growth rate modes interval (Sec. IV A). Whereas the time scale

*T*is defined as the period in which the perturbation growth is positive (i.e.,

_{WI}*ω*> 0). The variation of this non-dimensional average growth rate $ \sigma i \u0304 T W I $ as function of

_{i}*Re*

_{δ}is presented in Fig. 16. Although the underlying base flows are quite different (instability of a surface wave’s BBL vs. instability of internal waves in mid-water), the

*Re*

_{δ}= 350 limit of applicability of the quasi-steady linear stability analysis defined in this study occurs at a value of $ \omega i \u0304 T W I $ remarkably close to that computed by the studies mentioned above.

^{34,35}This agreement between the three different studies suggest that the particular instability criterion may be universal for a wider range of transient base flows.

## VI. CONCLUSIONS

In this study, a temporal instability map of the BBL under a solitary wave has been constructed. In this map, the connections between experimental observations, classical stability analysis, and fully non-linear numerical simulations have been established. The BBL under solitary wave was approximated by a flow in an oscillating water tunnel driven by a soliton-like wave induced pressure gradient similar to Sumer *et al.*^{9} Both quasi-steady/transient linear stability analysis and fully nonlinear two-dimensional simulations are carried out using high-order accuracy numerical methods.

The quasi-steady linear stability analysis has been employed to give an insight into the growth rate of the individual modes for each velocity profile according to a “momentary” criterion of instability as proposed Blondeaux *et al.*^{23} The *Re*_{δ} threshold above which the quasi-steady assumption holds is when the perturbation recovers its initial amplitude within the wave event. To formally identify this threshold, a transient linear stability analysis approach is adopted, which is also verified using fully non-linear DNS. This lower *Re*_{δ} limit corresponds to the case where the base flow is first destabilized before the flow reversal takes place.

The average perturbation growth rate corresponding to the above applicability threshold is then compared to the average growth rate associated with a similar instability criterion reported in the literature, albeit for different base flows, i.e., those driven in mid-water for progressive and solitary internal waves.^{34,35} Good agreement is found across all studies considered which suggests a possible universality of the particular instability criterion for a wider range of transient base flows.

After identifying the lower *Re*_{δ} limit of applicability of quasi-steady linear stability analysis, two-dimensional fully non-linear DNS have been used to classify the base flow into three possible regimes: unconditionally stable or laminar, conditionally unstable, and unconditionally unstable. Crucial to this classification is an alternative definition of “stability” which does not rely on perturbation growth rates. Instead, it is more qualitatively/visually based according to the laboratory observations of SU10. The map of instability illustrating the possible classifications of the base flow is shown in Fig. 17.

What is regarded as “stability” in the alternative, more visually oriented, flow classification approach is the ability of the perturbation to sufficiently grow in amplitude such that it can alter the base flow. In other words, the perturbation amplitude will be amplified, through linear instability, to the same order of magnitude of the base flow within the time scale of the wave. Vortex tubes then emerge as a result of the non-linear interaction of this amplified most unstable mode with the base flow. The critical *Re*_{δ} for the appearance of the vortex tubes, i.e., lower limit of the conditionally unstable regime, is investigated in more detail and the discrepancy with the experimental proposed limit is discussed. For a specific *Re*_{δ} case in the conditionally unstable regime, a very good agreement was found between numerical simulations and the SU10 measurements in terms of the bed shear stress and the characteristics, i.e., spacing and time of formation, of the formed vortex tubes. This lower limit of the unconditionally unstable regime is found to match the *Re*_{δ} limit for the transition to turbulence identified by Ozdemir *et al.*^{14}

The variation of the shape of the neutral curves with increasing *Re*_{δ} is also addressed in this study. An immediate consequence of such a variation is the existence of two primary modes of instability: post-$ M 1 $ and pre-$ M 2 $ flow reversal modes. The dominance of either of these modes is reflected in the resulting streamwise spacing between the generated vortex tubes. The post-flow reversal mode, which is related to the instability of the developed shear layer, is characterized by the highest growth rate and dominates for relatively low *Re*_{δ} flows, i.e., in the conditionally unstable regime. As the *Re*_{δ} is increased and the flow enters the second sub-regime of the unconditionally unstable regime, the *M*2 mode, which is longer in wavelength, destabilizes first and can grow sufficiently fast to dominate, leading to another type of instability different than the traditional shear layer instability. As *Re*_{δ} is increased further, *M*2 establishes a distinct signature within the neutral curve which starts to protrude into the acceleration phase. No inflection point exists within the base flow velocity profiles during this stage in flow evolution, which suggests a viscous origin of the instability linked to *M*2.

At the lower end of the *Re*_{δ} range where *M*2 is dominant when low initial amplitude perturbation perturbations are inserted at *ωt* = 0^{o}, a competition between the two modes takes place if one pushes further in time the above time of perturbation injection. If perturbation insertion is sufficiently delayed, the high growth of *M*1 can overcome the early destabilization of *M*2 and it is what sets the spacing between the generated vortex tubes. Effectively, another layer of base flow classification is therefore introduced, as shown in Fig. 17. For sufficiently high *Re*_{δ}, *M*2 always dominates regardless of perturbation insertion time. Finally, for the unconditionally unstable regime, increase of the amplitude of the insertion perturbation can lead to a non-linear transition to turbulence well before the flow reversal, a transition which cannot be captured using the two-dimensional approach employed in this study and which lies outside the scope of the current work.

This study has focused on the two-dimensional investigation of the stability of the BBL flow under solitary waves using both linear and non-linear approaches. The previously proposed, experimentally based, classification of the base flow as function of its Reynolds number has been revisited in a more detailed manner combined with an extension to larger values of this governing parameter. This study focuses on establishing a connection between the linear/non-linear analysis and the experimental observation from the point of view of the classical transition (two to three-dimensional) transition to turbulence, as observed in most free shear flows.^{15} Future work will focus on extending this study to account for the other types of transition^{37} that can take place in the place of this classical transition, namely, in a three-dimensional flow field. Moreover, a possible correlation between the wavelength of the coherent structure, linked to the dominant mode of either *M*1 or *M*2, and the resulting streamwise turbulence length scales. To this end, the findings of this study will be the foundation for a more detailed future three-dimensional stability analysis and direct numerical simulations. Additionally, the similarities between the results of this study and the near-bed instability under the internal solitary wave of elevation,^{38,39} the internal analog of a surface solitary wave, can be an additional area of future investigation. A topic of longer-term study is the implication of the developed instability and eventually transition to fully turbulent BL, on the near-bed sediment suspension and morphology.^{40}

## Acknowledgments

This work was supported by National Science Foundation Grant No. CBET-0756327 and CAREER award Grant (No. OCE-084558) and Bi-National U.S.-Israel Science Foundation Grant No. 2008051. The authors thank Professor M. Sumer for providing them with the experimental measurements. Two anonymous reviewers are thanked for a series of insightful comments.

### APPENDIX: ANALYTICAL SOLUTION FOR VISCOUS BOUNDARY LAYER FLOWS UNDER WEAKLY NONLINEAR TRANSIENT WAVE

This section summarizes the analytical solution of Liu and Orfilla^{5} for the viscous boundary layer flow under weakly nonlinear transient wave at a fixed location, i.e., *x* = 0. They showed that the rotational component of the boundary layer velocity can be written, in the dimension form, as follows:

where $ u w \u2217 $ denotes the free-stream velocity at the outer edge of the boundary layer (in this study, $ u w \u2217 t $ is simply $ U 0 \u2217 t \u2217 $ in (2)), *z*^{∗} is the vertical coordinate normal to the bed pointing upward, and *ν*^{∗} is the fluid kinematic viscosity. The corresponding laminar bed shear stress, $ \tau b \u2217 $, can be calculated as follows:

where *ρ* is the fluid density.

## REFERENCES

*Re*

_{δ}= 1100.