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.

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 Liu4 showed that the leading tsunami wave can be adequately modeled by a linear superposition of a number of sech2-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 Orfilla5 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 Cowen6,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(104) and the boundary layer flow was always laminar. Specifically, the free-stream Reynolds number is defined as

Re = a U o m ν ,
(1)

where U o m 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 , to the wave angular frequency ω (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 × 106 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 < 2 × 1 0 5 to laminar with two-dimensional vortex tubes near the bed 2 × 1 0 5 Re < 5 × 1 0 5 , and finally, to a transitional regime characterized by the appearance of single/multiple turbulent spots Re > 5 × 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 Orfilla5 and the previous wave tank measurements of Liu et al.6 for the laminar Re range (Re < 2 × 105). In the transitional regime, Re > 5 × 105, 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 × 106, 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 Blondeaux11 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 Blondeaux12 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 , 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 × 105. Moreover, the transition to turbulence was suggested to take place immediately above this cutoff, i.e., Re > 5 × 105.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 Scandura13 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 × 105). 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 Pedersen24 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.

Consider a two-dimensional solitary wave of height H that propagates in a constant water depth h. Hereafter, the sign indicates dimensional quantity. The Cartesian coordinate system x , z 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 / h 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:

U 0 ( t ) = U 0 m sech 2 ω t ,
(2)

where U o m is the wave amplitude defined as

U 0 m = g h H h
(3)

and g is the gravitational acceleration. The wave frequency ω is defined as

ω = 3 g H 4 h 2
(4)

and, hence, a wave period (T) can be defined as

T = 2 π ω ,
(5)

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

Hereafter, we will omit the from the product ω t 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 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

1 ρ p 0 x = U 0 t t = 2 U 0 m ω sech 2 ( ω t ) tanh ( ω t ) ,
(6)

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 and the corresponding pressure gradient normalized by the absolute maximum pressure gradient under a solitary wave.

FIG. 1.

Sketches of temporal variation of the normalized pressure gradient (dashed line) and free stream velocity (solid line) under a solitary wave at a fixed spatial location.

FIG. 1.

Sketches of temporal variation of the normalized pressure gradient (dashed line) and free stream velocity (solid line) under a solitary wave at a fixed spatial location.

Close modal

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

u t = ( u ) u 1 ρ p + ν 2 u ,
(7)
. u = 0 ,
(8)

where u 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 = ( x , z ) .

Consequently, the total pressure (p) in Eq. (7) can be divided into two parts,

p = p 0 x , t + p x , z , t ,
(9)

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 p 0 / x given in Eq. (6) that is constant in space yet variable in time.

To nondimensionalize the problem, we use U o m as a velocity scale and the Stokes laminar boundary layer thickness δ,

δ = 2 ν ω 1 2 ,
(10)

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

u = 0 ,
(11)
u t + ( u ) u = p + 1 Re δ 2 u ,
(12)

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

Re δ = δ U o m ν .
(13)

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

Re δ = 2 Re .
(14)

From Eqs. (3), (10), and (13), Reδ can be expressed as

Re δ = 16 3 H h h δ .
(15)

From the equation above, in contrast to previous studies,23 it is clear that the parameter space of the problem in hand, following the assumptions in (II A), collapses into a single parameter which is the Reynolds number Reδ.

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 Henningson26 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:

p = P ( x ) + p ̃ ( x , z , t ) ,
(16)
u = ( U ( z ) , 0 ) + ( u ̃ , w ̃ ) .
(17)

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:

u ̃ x + w ̃ z = 0 ,
(18)
u ̃ t + U u ̃ x + w U = p ̃ x + 1 Re δ 2 u ̃ ,
(19)
w ̃ t + U w ̃ x = p ̃ z + 1 Re δ 2 w ̃ .
(20)

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

t + U x 2 U x 1 Re δ 4 w ̃ = 0 .
(21)

The boundary conditions w ̃ = w ̃ / z = 0 both at the bed z = 0 and in the far-field z 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:

w ̃ = w ¯ ( z ) e i α x ω p t ,
(22)

where α is the wavenumber in the x direction ( α = 2 π λ x , λ x is the stream-wise wavelength of the perturbation, which eventually determines the spacing between the resulting vortex tubes), w ¯ is the perturbation wave-like amplitude, and ωp is the perturbation frequency ω p = ω r + i ω 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:

i ω p k 2 D z 2 w ¯ + L O S w ¯ = 0 ,
(23)

with the OS operator LOS defined as

L O S = i α U k 2 D z 2 + i α U + 1 Re δ k 2 D z 2 2 .
(24)

In the equation above, Dz is the derivative operator, with respect to z. k = α for a two-dimensional case. For the temporal stability analysis, ωp is complex and it appears as the eigenvalue in the OS equation for each profile. If ωi > 0, this indicates an unstable mode (i.e., growth in time), whereas ωi < 0 corresponds to a decaying mode.

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, ω t critical which has its corresponding critical wavenumber α critical . The maximum growth rate, ωimax, occurs at ω t ω i max and for a corresponding wavenumber α ω 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

w ̃ = w ̂ ( z , t ) e i α x .
(25)

Substituting Eq. (25) in Eq. (21), the general form of Eq. (23) is as follows:

k 2 D z 2 w ̂ t + L O S w ̂ = 0 ,
(26)

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

w ̂ t = L 1 n w ̂ ,
(27)

where L 1 = L O S / k 2 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 ̂ ( z , t ) over the domain volume V as follows:

q ( t ) = 1 2 V w ̂ 2 ( z , t ) d V .
(28)

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 3 .

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.

For the numerical solution of either the classical eigenvalue problem, Eq. (23), or the initial value problem, Eq. (26), a Chebyshev collocation method27 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 Orszag28 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),

i ω p w ¯ = k 2 D z 2 1 L O S w ¯ = L 1 w ¯ .
(29)

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 Δ ω t = 0 . 1 within which it is safe to assume that the base flow is nearly steady. The length of this interval, i.e., Δ ω t , can be translated into an equivalent perturbation growth time, i.e., tp, using the problem time scale i.e., δ / U 0 m . Within this Δ ω t interval, the linear operator L1 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 Δtp to proceed from step n to step n + 1 for the each of the total interval time tp as follows:

w ̂ n + 1 w ̂ n Δ t p = 1 2 L 1 n w ̂ n + 1 + L 1 n w ̂ n .

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.

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 Nx 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).

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.

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 × 104) 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 ωtc ≈ 40. The critical wave phase happens after the flow reversal which takes place at ωt ≈ 30 for all cases.

FIG. 2.

Growth contour of ωi > 0 for Reδ = 211 throughout the wave deceleration (thick solid line ωi = 0, contour value spacing Δωi = 0.001, and ωimax = 1.28 × 10−2). Horizontal dashed lines bound the modes with highest growth rates, i.e., α 0 . 43 , 0 . 47 . The table below the figure summarizes the following quantities: critical wave phase, critical wavenumber, maximum growth rate, wave phase of the maximum growth rate, and the wavenumber with the maximum growth rate (see also Sec. II B 2 for definitions).

FIG. 2.

Growth contour of ωi > 0 for Reδ = 211 throughout the wave deceleration (thick solid line ωi = 0, contour value spacing Δωi = 0.001, and ωimax = 1.28 × 10−2). Horizontal dashed lines bound the modes with highest growth rates, i.e., α 0 . 43 , 0 . 47 . The table below the figure summarizes the following quantities: critical wave phase, critical wavenumber, maximum growth rate, wave phase of the maximum growth rate, and the wavenumber with the maximum growth rate (see also Sec. II B 2 for definitions).

Close modal

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, α critical = 0 . 44 , as identified by the quasi-steady analysis. Across all Reδ, the highest growth rate modes, i.e., α ω i max , always lies within the range of α 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, q0, 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 ω t 4 , 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.

FIG. 3.

Time series for the normalized kinetic energy throughout the deceleration phase for Reδ = 211 for three wavenumbers: α = 0 . 43 dashed line , 0 . 45 solid line , and 0 . 47 dotted line .

FIG. 3.

Time series for the normalized kinetic energy throughout the deceleration phase for Reδ = 211 for three wavenumbers: α = 0 . 43 dashed line , 0 . 45 solid line , and 0 . 47 dotted line .

Close modal

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 ω t c r i t i c a l 3 0 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.

FIG. 4.

Times series for the normalized vertical perturbation kinetic energy for α = 0.45 throughout the deceleration phase for Reδ = 330 (dotted line), 350 (solid line), and 375 (dashed line).

FIG. 4.

Times series for the normalized vertical perturbation kinetic energy for α = 0.45 throughout the deceleration phase for Reδ = 330 (dotted line), 350 (solid line), and 375 (dashed line).

Close modal
FIG. 5.

Growth contour of ωi > 0 for Reδ = 350 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). Vertical dashed line indicates the start of the flow reversal. Table quantities are as defined in the caption of Fig. 2.

FIG. 5.

Growth contour of ωi > 0 for Reδ = 350 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). Vertical dashed line indicates the start of the flow reversal. Table quantities are as defined in the caption of Fig. 2.

Close modal

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.

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).

In Secs. IV B 1 and IV B 2, the results of the two-dimensional DNS analysis will be presented as an additional confirmation of this qualitatively/visually based classification which is consistent with SU10’s observations.

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 Blondeaux12 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 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 × 105) 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 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 

FIG. 6.

Growth contour of ωi > 0 for Reδ = 920 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). Table quantities are as defined in the caption of Fig. 2.

FIG. 6.

Growth contour of ωi > 0 for Reδ = 920 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). Table quantities are as defined in the caption of Fig. 2.

Close modal
FIG. 7.

Time series of the normalized bed shear stress for Reδ = 920 at two different spatial locations. Solid lines represent DNS results, whereas dashed lines represent laboratory measurements of SU10.

FIG. 7.

Time series of the normalized bed shear stress for Reδ = 920 at two different spatial locations. Solid lines represent DNS results, whereas dashed lines represent laboratory measurements of SU10.

Close modal

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 × 1 0 5 Re δ = 630 . In contrast, Vittori and Blondeaux11,12 identified a slightly higher Reynolds number limit, Re = 5 × 1 0 5 Re δ = 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 U0m).

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 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 × 1 0 4 Re δ = 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.

FIG. 8.

Time series of the normalized bed shear stress for Reδ = 450 for two case: the unperturbed case (i.e., analytical laminar solution, solid line) and perturbed case (dashed line).

FIG. 8.

Time series of the normalized bed shear stress for Reδ = 450 for two case: the unperturbed case (i.e., analytical laminar solution, solid line) and perturbed case (dashed line).

Close modal

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 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 12 to O 1 0 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 GW of each of the most unstable modes for the different Reδ cases is shown in Fig. 9. GW is the equivalent of the amplification factor 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 GW for dominant perturbation wavelength, i.e., λx ≈ 15, in the first sub-regime, i.e., 1900 < Reδ < 5000. The solid curves are representative of the transient variation of GW in the second sub-regime, i.e., Reδ > 5000, where the wavelength of the most unstable mode is λx ≈ 21. As shown in this figure, GW 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 Reδ increases, which agrees with the quasi-steady linear stability analysis (e.g., Figs. 5 and 6).

FIG. 9.

Transient growth as characterized by sample curves of amplification factor GW of the most unstable mode for the two categories of the unconditionally unstable regime. Dashed lines correspond to the growth of the λx ≈ 15 mode, whereas solid lines represent the λx ≈ 21 mode.

FIG. 9.

Transient growth as characterized by sample curves of amplification factor GW of the most unstable mode for the two categories of the unconditionally unstable regime. Dashed lines correspond to the growth of the λx ≈ 15 mode, whereas solid lines represent the λx ≈ 21 mode.

Close modal

As previously discussed, depending on Reδ, there are two distinct dominant instability modes. The first mode, hereafter regarded as M1, is associated with the highest energy growth rate, which consistently takes place after the flow reversal. The wavenumber of this mode is roughly around α 0 . 47 , 0 . 42 which corresponds to a non-dimensional wavelength of λ x 13 . 5 , 15 . Above a certain Reδ threshold, another mode, hereafter regarded M2, 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 M2 is much smaller than the highest energy growth rate associated with M1 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 M1 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 M1 observed in the two-dimensional DNS (e.g., Fig. 9). This numerical observation of M1 dominance fits well with the average distance between the successive vortex tubes observed in the SU10 study, i.e., λ Sumer 13 , 14 . 5 , and in previous numerical simulations12–14 within this particular Reδ range. In the first sub-regime of the unconditionally unstable regime, i.e., 1900 < Reδ < 5000, although M2 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 M1 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.

FIG. 10.

Variation of the neutral curves ω i = 0 as a function of Reδ. The solid ellipse represents the locus of centers of the maximum ωi of the different cases. The two vertical dashed lines correspond to the wave crest (left) and the flow reversal phase (right).

FIG. 10.

Variation of the neutral curves ω i = 0 as a function of Reδ. The solid ellipse represents the locus of centers of the maximum ωi of the different cases. The two vertical dashed lines correspond to the wave crest (left) and the flow reversal phase (right).

Close modal

As Reδ is increased further, the neutral curves start to show a visible protrusion towards a smaller wavenumber that of M2. Eventually, this protrusion penetrates into the acceleration phase region. As discussed earlier, for Reδ > 5000, M2 starts to grow at an earlier phase and it can actually overtake the high growth rate of M1 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 GW for the two modes M1 and M2. For this latter Reδ value case, M2 has the highest growth rate prior to the flow reversal, i.e., ωt < 30. After the flow reversal takes place, M1 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 M1 (i.e., λx ≈ 15; see also the left panel of Fig. 15 and its discussion).

FIG. 11.

Growth contour of ωi > 0 for Reδ = 4000 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). Table quantities are as defined in the caption of Fig. 2.

FIG. 11.

Growth contour of ωi > 0 for Reδ = 4000 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). Table quantities are as defined in the caption of Fig. 2.

Close modal
FIG. 12.

Transient growth of GW for Reδ = 4000. Dashed lines are for M1, i.e., λx ≈ 15, whereas solid lines are for M2, i.e., λx ≈ 21.

FIG. 12.

Transient growth of GW for Reδ = 4000. Dashed lines are for M1, i.e., λx ≈ 15, whereas solid lines are for M2, i.e., λx ≈ 21.

Close modal
FIG. 15.

Visualization of the perturbation vorticity field (normalized by the corresponding maximum value) for Reδ = 4000 (Figs. 15(a)–15(c) and Reδ = 8000 (Figs. 15(d)–15(f)) at different wave phases ωt. First row: linear stability analysis. Second and third rows: DNS results. In the DNS results, only a portion of the computational domain is shown. In the left column, M1 is the dominant instability mode. In the right column, M2 dominates. The horizontal and vertical coordinates in each panel are normalized by the boundary layer thickness δ.

FIG. 15.

Visualization of the perturbation vorticity field (normalized by the corresponding maximum value) for Reδ = 4000 (Figs. 15(a)–15(c) and Reδ = 8000 (Figs. 15(d)–15(f)) at different wave phases ωt. First row: linear stability analysis. Second and third rows: DNS results. In the DNS results, only a portion of the computational domain is shown. In the left column, M1 is the dominant instability mode. In the right column, M2 dominates. The horizontal and vertical coordinates in each panel are normalized by the boundary layer thickness δ.

Close modal

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: ωtins = 0,  18,  and 30 are also included in this figure.

FIG. 13.

Growth contour of ωi > 0 for Reδ = 6300 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). The three vertical dashed lines mark the different insertion phases. Table quantities are as defined in the caption of Fig. 2.

FIG. 13.

Growth contour of ωi > 0 for Reδ = 6300 throughout the wave deceleration (thick solid line ωi = 0 and contour value spacing Δωi = 0.001). The three vertical dashed lines mark the different insertion phases. Table quantities are as defined in the caption of Fig. 2.

Close modal

At the same value of Reδ, Fig. 14 shows the transient variation of GW for each M1 and M2 for the three different times of perturbation insertion i.e. , ω t ins = 0 , 1 8 , and 3 0 where the initial relative amplitude of the inserted perturbation is O 1 0 12 . In the calculation of GW for the different cases of insertion, the kinetic energy of the vertical velocity component for the particular mode (i.e., q W ω t ) is normalized by its initial value at the insertion time q W ω t ins . In the left panel of Fig. 14, M2 starts to grow immediately at ωt = 0, as the adverse pressure gradient is established, with a relatively high growth rate compared to M1. M1 picks up around the flow reversal phase, i.e., ωt ≈ 30, but cannot overcome the earlier growth of M2 (as in the case Reδ = 4000, Fig. 12) provided that the perturbation was inserted prior to flow reversal i.e., ω t 3 0 . As the perturbation is inserted closer to the flow reversal phase angle, the most unstable mode changes from M2 to M1 as shown for the ωtins ≈ 30 case.

FIG. 14.

Transient growth of GW for Reδ = 6300 for three different noise insertion phase cases: ωtins = 0 (a), ωtins = 18 (b), and ωtins = 30 (c). M1 is represented by the dashed line, whereas solid lines correspond to M2.

FIG. 14.

Transient growth of GW for Reδ = 6300 for three different noise insertion phase cases: ωtins = 0 (a), ωtins = 18 (b), and ωtins = 30 (c). M1 is represented by the dashed line, whereas solid lines correspond to M2.

Close modal

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. M1 dominates for Reδ = 4000, whereas M2 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.

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 Koseff34 and Barad and Fringer35 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 σ i ̄ T w > 5 , where σ i ̄ is the growth rate of the instabilities averaged over a time interval Tw during which a parcel of fluid within stratified water is subjected to a wave-driven background vertical shear corresponding to Richardson number36 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 ω i ̄ T W I may be calculated where ω i ̄ is average positive growth of the critical wavenumber, i.e., α critical , found to always lie within the highest growth rate modes interval (Sec. IV A). Whereas the time scale TWI is defined as the period in which the perturbation growth is positive (i.e., ωi > 0). The variation of this non-dimensional average growth rate σ i ̄ T W I as function of 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 ω i ̄ 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.

FIG. 16.

Non-dimensional average growth rate variation as a function of Reδ. Stable, ● and unstable, ■. Horizontal dashed line represents the applicability limit of the quasi-steady linear stability analysis.

FIG. 16.

Non-dimensional average growth rate variation as a function of Reδ. Stable, ● and unstable, ■. Horizontal dashed line represents the applicability limit of the quasi-steady linear stability analysis.

Close modal

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.

FIG. 17.

Instability map of the BBL flow under solitary waves. The flow is classified in accordance to four criteria: quasi-steady steady linear stability analysis (upper bar), transient linear stability analysis (second bar from top), appearance of the vortex tubes (third from top), and dominant mode of instability (bottom bar).

FIG. 17.

Instability map of the BBL flow under solitary waves. The flow is classified in accordance to four criteria: quasi-steady steady linear stability analysis (upper bar), transient linear stability analysis (second bar from top), appearance of the vortex tubes (third from top), and dominant mode of instability (bottom bar).

Close modal

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 M2 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, M2 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 M2.

At the lower end of the Reδ range where M2 is dominant when low initial amplitude perturbation perturbations are inserted at ωt = 0o, 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 M1 can overcome the early destabilization of M2 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δ, M2 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 transition37 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 M1 or M2, 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 

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.

This section summarizes the analytical solution of Liu and Orfilla5 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:

u z , t = z 4 π 0 t u w t t τ 3 exp z 2 4 ν t τ d τ ,
(A1)

where u w denotes the free-stream velocity at the outer edge of the boundary layer (in this study, u w t is simply U 0 t 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, τ b , can be calculated as follows:

τ b t ρ = ν π 0 t u w t / τ t τ 3 d τ ,
(A2)

where ρ is the fluid density.

1.
P. L.-F.
Liu
,
P.
Lynett
,
H.
Fernando
,
B. E.
Jaffe
,
H.
Fritz
,
B.
Higman
,
R.
Morton
,
J.
Goff
, and
C.
Synolakis
, “
Observations by the international tsunami survey team in Sri Lanka
,”
Science
308
,
1595
(
2005
).
2.
H.
Yeh
,
M.
Francis
,
C.
Peterson
,
T.
Katada
,
G.
Latha
,
R.
Chadha
,
J.
Singh
, and
G.
Raghuraman
, “
Effects of the 2004 great sumatra tsunami: Southeast Indian coast
,”
J. Waterw., Port, Coastal, Ocean Eng.
133
,
382
(
2007
).
3.
P. A.
Madsen
,
D. R.
Fuhrman
, and
H. A.
Schäffer
, “
On the solitary wave paradigm for tsunamis
,”
J. Geophys. Res.: Oceans
113
,
C12012
, doi:10.1029/2008jc004932 (
2008
).
4.
I.
Chan
and
P. L.-F.
Liu
, “
On the runup of long waves on a plane beach
,”
J. Geophys. Res.: Oceans
117
,
C08006
, doi:10.1029/2012jc007994 (
2012
).
5.
P. L.-F.
Liu
and
A.
Orfila
, “
Viscous effect on transient long-wave propagation
,”
J. Fluid Mech.
520
,
83
(
2004
).
6.
P. L.-F.
Liu
,
Y. S.
Park
, and
E. A.
Cowen
, “
Boundary layer flow and bed shear stress under a solitary wave
,”
J. Fluid Mech.
574
,
449
(
2007
).
7.
Y. S.
Park
,
J.
Verschaeve
,
G. K.
Pedersen
, and
P. L.-F.
Liu
, “
Boundary-layer flow and bed shear stress under a solitary wave: Revision
,”
J. Fluid Mech.
753
,
554
(
2014
).
8.
W. O.
Criminale
,
T. L.
Jackson
, and
R. D.
Joslin
,
Theory and Computation of Hydrodynamic Stability
(
Cambridge University Press
,
2003
).
9.
B. M.
Sumer
,
P. M.
Jensen
,
L. B.
Sørensen
,
J.
Fredsøe
,
P. L.-F.
Liu
, and
S.
Carstensen
, “
Coherent structures in wave boundary layers. Part 2. Solitary motion
,”
J. Fluid Mech.
646
,
207
(
2010
); Referred to in text as SU10.
10.
S.
Carstensen
,
B. M.
Sumer
, and
J.
Fredsøe
, “
Coherent structures in wave boundary layers. Part 1. Oscillatory motion
,”
J. Fluid Mech.
646
,
169
(
2010
).
11.
G.
Vittori
and
P.
Blondeaux
, “
Turbulent boundary layer under a solitary wave
,”
J. Fluid Mech.
615
,
433
(
2008
).
12.
G.
Vittori
and
P.
Blondeaux
, “
Characteristics of the boundary layer at the bottom of a solitary wave
,”
Coastal Eng.
58
,
206
(
2011
).
13.
P.
Scandura
, “
Two-dimensional vortex structures in the bottom boundary layer of progressive and solitary waves
,”
J. Fluid Mech.
728
,
340
(
2013
).
14.
C. E.
Ozdemir
,
T.-J.
Hsu
, and
S.
Balachandar
, “
Direct numerical simulations of instability and boundary layer turbulence under a solitary wave
,”
J. Fluid Mech.
731
,
545
(
2013
).
15.
P. G.
Drazin
and
W. H.
Reid
,
Hydrodynamic Stability
(
Cambridge University Press
,
2004
).
16.
T.
Herbert
, “
Parabolized stability equations
,”
Annu. Rev. Fluid Mech.
29
,
245
(
1997
).
17.
P. J.
Blennerhassett
and
A. P.
Bassom
, “
The linear stability of flat stokes layers
,”
J. Fluid Mech.
464
,
393
(
2002
).
18.
F. J.
Poulin
,
G. R.
Flierl
, and
J.
Pedlosky
, “
Parametric instability in oscillatory shear flows
,”
J. Fluid Mech.
481
,
329
(
2003
).
19.
J.
Luo
and
X.
Wu
, “
On the linear instability of a finite stokes layer: Instantaneous versus floquet modes
,”
Phys. Fluids
22
,
054106
(
2010
).
20.
S.
Shen
, “
Some considerations on the laminar stability of time-dependent basic flows
,”
AIAA J.
28
,
397
(
1961
).
21.
G.
Seminara
and
P.
Hall
, “
Linear stability of slowly varying unsteady flows in a curved channel
,”
Proc. R. Soc. London, Ser. A
346
,
279
(
1975
).
22.
P.
Hall
and
K.
Parker
, “
The stability of the decaying flow in a suddenly blocked channel
,”
J. Fluid Mech.
75
,
305
(
1976
).
23.
P.
Blondeaux
,
J.
Pralits
, and
G.
Vittori
, “
Transition to turbulence at the bottom of a solitary wave
,”
J. Fluid Mech.
709
,
396
(
2012
).
24.
J. C.
Verschaeve
and
G. K.
Pedersen
, “
Linear stability of boundary layers under solitary waves
,”
J. Fluid Mech.
761
,
62
(
2014
).
25.
R.
Grimshaw
, “
The solitary wave in water of variable depth. Part 2
,”
J. Fluid Mech.
46
,
611
(
1971
).
26.
P. J.
Schmid
and
D. S.
Henningson
,
Stability and Transition in Shear Flows
(
Springer
,
2001
), Vol.
142
.
27.
J. P.
Boyd
,
Chebyshev and Fourier Spectral Methods
(
Courier Dover Publications
,
2013
).
28.
S. A.
Orszag
, “
Accurate solution of the orr–sommerfeld stability equation
,”
J. Fluid Mech.
50
,
689
(
1971
).
29.
P. J.
Diamessis
,
J.
Domaradzki
, and
J.
Hesthaven
, “
A spectral multidomain penalty method model for the simulation of high Reynolds number localized incompressible stratified turbulence
,”
J. Comput. Phys.
202
,
298
(
2005
).
30.
P. J.
Diamessis
and
L. G.
Redekopp
, “
Numerical investigation of solitary internal wave-induced global instability in shallow water benthic boundary layers
,”
J. Phys. Oceanogr.
36
,
784
(
2006
).
31.
G. E.
Karniadakis
,
M.
Israeli
, and
S. A.
Orszag
, “
High-order splitting methods for the incompressible Navier-Stokes equations
,”
J. Comput. Phys.
97
,
414
(
1991
).
32.
D.
Gottlieb
and
J. S.
Hesthaven
, “
Spectral methods for hyperbolic problems
,”
J. Comput. Appl. Math.
128
,
83
(
2001
).
33.
See supplementary material at http://dx.doi.org/10.1063/1.4916560 for movies showing the vorticity evolution forReδ = 1100.
34.
C. D.
Troy
and
J. R.
Koseff
, “
The instability and breaking of long internal waves
,”
J. Fluid Mech.
543
,
107
(
2005
).
35.
M. F.
Barad
and
O. B.
Fringer
, “
Simulations of shear instabilities in interfacial gravity waves
,”
J. Fluid Mech.
644
,
61
(
2010
).
36.
S. A.
Thorpe
,
The Turbulent Ocean
(
Cambridge University Press
,
2005
).
37.
P. J.
Schmid
, “
Nonmodal stability theory
,”
Annu. Rev. Fluid Mech.
39
,
129
(
2007
).
38.
M.
Stastna
and
K. G.
Lamb
, “
Vortex shedding and sediment resuspension associated with the interaction of an internal solitary wave and the bottom boundary layer
,”
Geophys. Res. Lett.
29
,
7
, doi:10.1029/2001GL014070 (
2002
).
39.
M.
Stastna
and
K. G.
Lamb
, “
Sediment resuspension mechanisms associated with internal waves in coastal waters
,”
J. Geophys. Res.: Oceans
113
,
C10016
, doi:10.1029/2007jc004711 (
2008
).
40.
Y.-J.
Chou
and
O. B.
Fringer
, “
A model for the simulation of coupled flow-bed form evolution in turbulent flows
,”
J. Geophys. Res.: Oceans
115
,
C10041
, doi:10.1029/2010jc006103 (
2010
).

Supplementary Material