Roughness-induced boundary-layer instability beneath internal solitary waves

Abstract

We investigate the effects of bottom roughness on bottom boundary-layer (BBL) instability beneath internal solitary waves (ISWs) of depression. Applying both two-dimensional (2-D) numerical simulations and linear stability theory, an extensive parametric study explores the effect of the Reynolds number, pressure gradient, roughness (periodic bump) height $h_b$ and roughness wavelength $\lambda _b$ on BBL instability. The simulations show that small $h_b$, comparable to that of laboratory-flume materials ($\sim$100 times less than the thickness of the viscous sublayer $\delta _v$), can destabilize the BBL and trigger vortex shedding at critical Reynolds numbers much lower than what occurs for numerically smooth surfaces. We identify two mechanisms of vortex shedding, depending on $h_b/\delta _v$. For $h_b/\delta _v \gtrapprox 1$, vortices are forced directly by local flow separation in the lee of each bump. Conversely, for $h_b/\delta _v \lessapprox 10^{-1}$ the roughness seeds perturbations in the BBL, which are amplified by the BBL flow. Roughness wavelengths close to those associated with the most unstable BBL mode, as predicted by linear instability theory, are preferentially amplified. This resonant amplification nature of the BBL flow, beneath ISWs, is consistent with what occurs in a BBL driven by surface solitary waves and by periodic monochromatic waves. Using the $N$-factor method for Tollmien–Schlichting waves, we propose an analogy between the roughness height and seed noise required to trigger instability. Including surface roughness, or more generally an appropriate level of seed noise, reconciles the discrepancies between the vortex-shedding threshold observed in the laboratory versus that predicted by otherwise smooth-bottomed 2-D spectral simulations.

1. Introduction

Internal solitary waves (ISWs) are nonlinear waves of large amplitude that can propagate tens to hundreds of kilometres along the pycnoclines of lakes and coastal oceans. As they shoal, ISWs of depression interact with the bottom boundary-layer (BBL) (Trowbridge & Lentz 2018) by imposing a streamwise velocity (pressure) distribution that is maximum (minimum) under the trough. The ISW-imposed pressure leads to rapid, unsteady acceleration and deceleration of the BBL flow beneath the front and rear shoulders of the ISW, respectively (e.g. Boegman & Stastna 2019; Zulberti et al. 2020). For large-amplitude ISWs propagating over a flat bottom, laboratory experiments (Carr et al. 2008; Aghsaee & Boegman 2015; Zahedi et al. 2021) and numerical simulations (Diamessis & Redekopp 2006; Aghsaee et al. 2012; Sakai et al. 2020; Ellevold & Grue 2023; Posada-Bedoya et al. 2024) have demonstrated that ISW-induced currents can lead to hydrodynamic instabilities in the BBL. At laboratory scales ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\sim \mathcal {O}(100)$ , defined in § 1.2), the instability is convective (Posada-Bedoya et al. 2024), as it falls behind the ISW. The instabilities in the BBL may be amplified by the flow (Huerre & Monkewitz 1990; Schmid & Henningson 2001), depending on the magnitude and energy distribution of the background seed noise (Verschaeve & Pedersen 2014; Posada-Bedoya et al. 2024) as well as the characteristics of the BBL flow (Reynolds number and pressure gradient) (Aghsaee et al. 2012).

Parameterizing and understanding the mechanisms leading to BBL instability beneath ISWs is crucial for closing ISW energy budgets (Moum et al. 2007b ; Zahedi et al. 2021) and predicting sediment transport processes (Boegman & Stastna 2019). Internal solitary waves carry significant energy shoreward (Chang et al. 2006; Shroyer et al. 2010), and as they approach the continental shelf, their energy is dissipated largely as a result of strong bottom interaction (Boegman et al. 2005; Moum et al. 2007a ; Aghsaee et al. 2010). Therefore, the stability of the BBL beneath ISWs will determine dissipation (Klymak & Moum 2003; Shroyer et al. 2009) and propagation distances (Shroyer et al. 2010; Zahedi et al. 2021, 2023): unstable waves have dissipative length scales of $\sim$ 100 wavelengths, whereas stable waves propagate significantly further, $\sim$ 1000 wavelengths. Predicting the BBL stability is also important for closing the turbulent kinetic energy budget in the BBL. Enhanced near-bed shear, induced by the passage of ISWs, has been shown to increase both production and dissipation rate in the BBL by three orders of magnitude (Zulberti et al. 2020, 2022). This, in turn, has direct implications for the sediment resuspension and transport (Stastna & Lamb 2008; Boegman & Stastna 2019).

Part of the difficulty in understanding the transition from laminar to turbulent flow lies in the number of factors that affect the transition, the most important ones being the pressure distribution in the external flow, the wall roughness and the characteristics of disturbances in the free-stream flow (Schlichting 1968). Previous studies on ISW-induced BBL instability have focused on the distribution of the external flow, by varying the pressure gradient and Reynolds number imposed by the ISW on the BBL (e.g. Diamessis & Redekopp 2006; Aghsaee et al. 2012; Ellevold & Grue 2023; Posada-Bedoya et al. 2024). In comparison, research on the roles of bottom roughness (e.g. Stastna & Lamb 2008; Carr et al. 2010) and of free-stream disturbances (e.g. Posada-Bedoya et al. 2024), in the ISW-induced BBL, remains scarce. Here, we focus on the role of bottom roughness in BBL instability beneath ISWs of depression.

1.1. Roughness-induced BBL instability

Simulations and laboratory experiments on the interaction of ISWs with rough bottom topography have considered topography with scales significantly larger than the induced boundary-layer thickness (Stastna & Lamb 2008; Carr et al. 2010; Olsthoorn & Stastna 2014; Harnanan et al. 2015,2017). For example, Stastna & Lamb (2008) simulated the propagation of ISWs of elevation over topography modelled as the superposition of three sinusoids with a maximum roughness height ( $h_b$ ) relative to the total depth $H$ , $h_b/H=0.03$ . Carr et al. (2010) conducted two-dimensional (2-D) simulations and laboratory experiments of ISWs propagating over sinusoidal corrugated beds. They report the smallest topography in the literature beneath ISWs of depression, with $h_b/H\approx 10^{-2}$ . In comparison, typical laboratory-flume bed materials for experiments on ISW-induced BBL instability including acrylic sediments (Aghsaee & Boegman 2015), smooth concrete (Zahedi et al. 2021) or glass (Carr et al. 2008), are expected to have relative roughness heights $h_b/H\sim 10^{-3}{-}10^{-6}$ . Based upon our simulations, these can be 10–100 times smaller than the laminar viscous sublayer thickness ( $\delta _v$ ) of the ISW-induced BBL. Despite the known presence of surface roughness in laboratory flumes, and its known theoretical relevance for BBL stability (Schlichting 1968), the physical mechanisms underlying such small-scale roughness-induced transition beneath ISWs remain unknown.

The roughness-induced BBL instability beneath ISWs investigated here is similar to the related problem of roughness-induced BBL instability beneath surface solitary waves (SSWs) (Sumer et al. 2010; Scandura 2013) and surface periodic waves (Blondeaux & Vittori 1994). Wall roughness has been shown to introduce perturbations to the surface-wave-driven BBL that are preferentially amplified at the wavelength of the most unstable mode of the BBL, as predicted by linear stability theory (Blondeaux & Vittori 1994; Scandura 2013). The similarity of the BBL flow beneath SSWs and ISWs raises the question of whether the BBL exhibits the same resonator behaviour beneath ISWs. Both types of waves (SSWs and ISWs) impose an unsteady, isolated pair of favourable and adverse pressure gradients on the BBL, with an inflection point in the velocity profile upstream of the BBL separation (Sumer et al. 2010; Aghsaee & Boegman 2015). Also, for both types of waves, Verschaeve & Pedersen (2014) predicted noise-amplifier behaviour of the BBL, with the stability determined by the accumulated amplification of perturbations of small initial amplitude (potentially introduced by small-scale bottom roughness). However, for both types of BBL flow, different studies reported different critical Reynolds numbers associated with instability. This is potentially explained by differences in the background noise level (Verschaeve & Pedersen 2014).

Figure 1. Stability diagrams in (a) $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space and (b) $a/H$ versus $Re_w$ space, showing our simulated ISW cases (numbered yellow and pink squares), the laboratory observations by Zahedi et al. (2021) (Z21) (circles) and Carr et al. (2008) (C08) (triangles) and the unstable ISW numerically simulated by Sakai et al. (2020) (S20) (blue square). Pink squares highlight the selected ISW cases to focus on the description of the results. The ISW numeration corresponds with the one in table 1. Black solid lines are stability curves from spectral 2-D simulations by (a) Aghsaee et al. (2012) (A12) and (b) Diamessis & Redekopp (2006) (D06). Blue and red lines are stability curves from 2-D simulations by Ellevold & Grue (2023) (E23) for two different pycnocline depths ( $d/H$ ). Blue and red markers are the associated laboratory experiments in Carr et al. (2008) selected by Ellevold & Grue (2023) to fit each stability curve. The C08, Z21 and E23 ISWs were generated by lock release, whereas the D06, A12, S20 and present study ISWs were generated by the solution of the Korteweg-de Vries (KdV) or DJL equations.

1.2. Vortex-shedding thresholds beneath ISWs

Using 2-D spectral simulations, Aghsaee et al. (2012) argued that BBL instability under an ISW of depression was determined by the non-dimensional pressure gradient ( $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ ) and the momentum thickness Reynolds number ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ ) at the separation point beneath the wave (figure 1), where

(1.1) \begin{equation} P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}= \left ( U_2+c \right ) \frac {U_2}{L_w{g}'}, \end{equation}

(1.2) \begin{equation} Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}= U_2 \sqrt {\frac {L_w}{\nu \left ( U_2+c \right )}}. \end{equation}

Here, $U_2$ is the absolute value of the maximum horizontal velocity at the wave trough, $c$ is the solitary wave phase speed, ${g}'=g\Delta \rho /\rho _0$ is the reduced gravity, $g$ is the gravitational constant, $\rho _0\,=1000\,\rm kg\,m^-{^3}$ is a reference density, $\Delta \rho$ is the density difference across the pycnocline, $\nu$ is the kinematic viscosity and $L_w$ is the horizontal wavelength scale (Michallet & Ivey 1999)

(1.3) \begin{equation} L_w=\frac {1}{a}\int _{-\infty }^{\infty }\eta _p(x){\rm d}x, \end{equation}

where $\eta _p(x)$ is the vertical displacement of the pycnocline and $a$ is the wave amplitude. More recently, Ellevold & Grue (2023) conducted 2-D simulations to reproduce the laboratory experiments by Carr et al. (2008). They proposed that, in addition to $a/H$ and $Re_w=c_0H/\nu$ , the pycnocline depth ( $d$ ) is a relevant parameter for BBL stability. They defined a critical threshold of the form $a/H=a_0(Re_w/Re_{w_{0}})^{-m_1}$ , where the parameters $a_0$ and $m_1$ depend on the relative depth of the pycnocline ( $d/H$ ). Here, $c_0$ is the linear wave phase speed.

At laboratory scales ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W} \sim \mathcal {O}(100)$ ), there is a discrepancy in the criteria for vortex shedding between experiments and spectral 2-D simulations (Aghsaee et al. 2012; Zahedi et al. 2021; Posada-Bedoya et al. 2024). Lock-release-generated ISWs generated in the laboratory (Carr & Davies 2006; Carr et al. 2008; Zahedi et al. 2021) and simulated by finite-volume 2-D solvers (Ellevold & Grue 2023) agree on the instability threshold, $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\approx 200$ (figure 1 a). However, 2-D spectral simulations of Dubreil–Jacotin–Long (DJL)-initialized ISWs (Aghsaee et al. 2012; Posada-Bedoya et al. 2024) predict a critical threshold at higher ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ ) (figure 1 a). It remains unknown why the 2-D spectral simulations have a different stability threshold from that observed in the laboratory.

The convective noise-amplifier nature of the BBL (Verschaeve & Pedersen 2014; Posada-Bedoya et al. 2024) supports the hypothesis that wall roughness might provide a mechanism for introducing seed perturbations to the BBL, which are convectively amplified beneath the ISW. Roughness occurs naturally in the laboratory but is not explicitly considered in idealized smooth-bed numerical simulations (e.g. Aghsaee et al. 2012; Sakai et al. 2020; Posada-Bedoya et al. 2024). Therefore, roughness, as a source of perturbations, is a potential feature that might explain the discrepancies between the experimental and spectral 2-D numerical results. To this point, Verschaeve & Pedersen (2014) explain why Sumer et al. (2010), Vittori & Blondeaux (2011) and Ozdemir et al. (2013) each obtained different critical Reynolds numbers for the instability of the BBL beneath SSWs due to differences in background noise across their numerical and experimental conditions. A similar explanation, with roughness as the distinguishing feature of differences in background noise, can be argued for the ISW case. Related to the present study, Scandura (2013) simulated SSW propagation over random bottom roughness with average height $h_b\approx 5\times 10^{-5}\,\rm m$ to reproduce the laboratory observations by Sumer et al. (2010).

1.3. Objectives

The objective of the present study is to investigate how bottom roughness may seed BBL instability beneath ISWs of depression, at scales comparable to the roughness of typical materials encountered in laboratory flumes. Using 2-D simulations and linear stability theory, we conducted an extensive parametric study on the effects of resolved sinusoidal roughness parameters (height and wavelength) on BBL instability beneath ISWs. We described the destabilization mechanisms and computed instability thresholds for different roughness heights. By interpreting BBL–roughness interaction as a mechanism for seeding perturbations within the BBL, we predicted the background noise levels associated with the different roughness heights. Finally, we discussed the relevance of our results for more realistic random roughness heights.

Figure 2. (a) Sketch of the domain, boundary conditions (b.c.), and the non-dimensional horizontal velocity field for the set-up of the numerical simulations. (b) Inset zoom of the bottom roughness elements indicates the roughness parameters ( $\lambda _b$ , $h_b$ ). (c) Inset zoom of the near-bed velocity indicates separation point ( $x_{sep}$ ), and boundary-layer parameters $\delta _v$ and $\delta _s$ . As $u_*=0$ at $x_{sep}$ , values of $5\nu /u_*$ are only shown upstream of $x_{sep}$ . Black lines show selected velocity profiles every $0.5x/H$ , with the zero of each profile indicated with a vertical dashed line. Profiles were scaled to fit the $0.5x/H$ spacing.

2. Methods

2.1. Problem definition

We conducted 2-D simulations of ISWs propagating over a flat bottom and a region of sinusoidal roughness. The geometry of the sinusoidal roughness was determined by two parameters: height $h_b$ and wavelength $\lambda _b$ (figure 2 b). The computational domain had a maximum depth $H$ and a length $L=16H$ . The ISWs were initialized by numerically solving the DJL equation using the algorithm of Turkington et al. (1991) as implemented by Dunphy et al. (2011). The initial condition was given by the DJL solution with the ISW trough at $x=0.25L$ . As in Carr et al. (2010), the present domain comprised a section of flat bottom, followed by rough topography, to allow for illustration of the effects of roughness on BBL instability. The region of sinusoidal roughness was located at a distance of $3L_w$ from the initial ISW trough, sufficient for the BBL beneath the initially inviscid DJL wave to develop, before interacting with the roughness region. We verified that the results were insensitive to further increasing this distance (not shown). The length of the roughness region was $0.31L$ , defined to accommodate 10 wavelengths of the largest simulated $\lambda _b$ (table 1) while balancing the computational demand with the total domain length. The results were insensitive to further increasing the length of the roughness region (not shown). A schematic of the problem is presented in figure 2(a, b).

Table 1. Parameters of simulated ISWs. In all cases, $h_{pyc}/H=0.02$ and $\Delta \rho /\rho _0=0.8$ . Boundary-layer thickness measures, $\delta _s$ and $\delta _v$ are scaled by the largest roughness height simulated $h_{b_0}/H=10^{-3}$ . For each ISW (i.e. each row, excluding cases A and B), four simulations were conducted with different bottom roughness (80 simulations): $h_b/H=10^{-3}$ , $h_b/H=10^{-4}$ , $h_b/H=10^{-5}$ and $h_b/H=0$ (flat bottom). ^†Cases A and B were simulated to investigate the effect of roughness wavelength and only considered $h_b/H=10^{-3}$ and $h_b/H=10^{-4}$ (16 simulations). A total of 96 simulations were conducted.

The roughness region was initiated and terminated with a minimum value of the sinusoidal function that defined the topography (see figure 2 b). This created a smooth transition between the roughness region and the flat bottom. Therefore, even though the bottom slope is a continuous function, the second-order derivative of the topography is discontinuous at the flat–rough bottom transition. Hence, the flow at the edges of the roughness region was different from the flow over the periodic bumps. As we are focused on the BBL interaction with the periodic topography, the flow near the flat–rough transitions was not considered here in detail.

The ISWs propagated along a quasi-two-layer density stratification defined via a hyperbolic tangent profile, widely used in numerical studies

(2.1) \begin{equation} \bar {\rho }(z)=\rho _0+\Delta \rho \tanh \left ( \frac {z-z_{pyc}}{h_{pyc}}\right ), \end{equation}

which represents a pycnocline of half-thickness $h_{pyc}$ , centred at depth $z_{pyc}$ , with $z$ positive upwards. In all cases, $h_{pyc}/H=0.02$ and $\Delta \rho /\rho _0=0.8$ .

Due to the streamwise variability of the BBL flow, boundary-layer parameters varied in the streamwise direction (figure 2 c). For scaling purposes, we defined the boundary-layer thickness representative of each ISW ( $\delta _{s}$ ) as the height of the 99% free-stream velocity ( $\delta _{99}$ ) at the ISW trough, where $U_{\infty }=U_2$ . The viscous sublayer thickness representative of each ISW ( $\delta _v$ ), was defined as the minimum value of the laminar sublayer height $5\frac {\nu }{u_{*}}$ upstream of the separation point, which is the region of the BBL where the flow was unstable. In figure 2(c), values of $5\nu /u_*$ are only shown upstream of $x_{sep}$ because $5\nu /u_* \to \infty$ as $x\to x_{sep}$ (because $u_*\to 0$ ). Here, $u_{*}$ is the friction velocity and $\tau _b$ is the wall shear stress, both varying in the streamwise direction and defined as

(2.2) \begin{equation} u_{*}(x)=\left (\frac {\tau _b}{\rho _0} \right )^{1/2}, \qquad \tau _b(x)=\nu \rho _0 {\left ( \frac {\partial u}{\partial z}\right )}_{z=0}. \end{equation}

2.2. Numerical model and simulation set-up

The 2-D numerical simulations were conducted with the pseudospectral code SPINS (Subich et al. 2013). SPINS solves the 2-D incompressible Navier–Stokes equations under the Boussinesq approximation

(2.3) \begin{gather} \frac {\partial {u}}{\partial {t}} + u\frac {\partial {u}}{\partial {x}} + w\frac {\partial {u}}{\partial {z}}=-\frac {1}{\rho _0}\frac {\partial {p}}{\partial {x}}+\nu \nabla ^2 u, \end{gather}

(2.4) \begin{gather} \frac {\partial {w}}{\partial {t}} + u\frac {\partial {w}}{\partial {x}} + w\frac {\partial {w}}{\partial {z}}=-\frac {1}{\rho _0}\frac {\partial {p}}{\partial {z}}+\nu \nabla ^2 w-\frac {\rho g}{\rho _0}, \end{gather}

(2.5) \begin{gather} \frac {\partial {\rho }}{\partial {t}} + u\frac {\partial {\rho }}{\partial {x}} + w\frac {\partial {\rho }}{\partial {z}}=\kappa \nabla ^2\rho , \end{gather}

(2.6) \begin{gather} \frac {\partial {u}}{\partial {x}} + \frac {\partial {w}}{\partial {z}} = 0 , \end{gather}

where $(x,z)$ are the horizontal and vertical coordinates, $(u,w)$ are the associated velocity vectors, $t$ is time, $p$ is the pressure, $\rho$ is the fluid density and $\kappa$ is the molecular diffusivity. For all cases, the simulated $Pr=\nu /\kappa =1$ .

SPINS implements a spectral collocation method for the simulation of the stratified Navier–Stokes equations in rectilinear and smooth curvilinear geometries. To solve for curvilinear geometries, the model defines a terrain-following grid in the physical space, which is mapped to a rectangular computation domain (see for more details Subich et al. 2013). The model solves the equations in the computational box and then maps the solution back to physical coordinates. The spectral accuracy in conjunction with a curvilinear geometry makes SPINS unique and suitable for the study of boundary-layer instability of stratified flows over non-flat-bottom boundaries. For the no-slip boundary condition (b.c.), the model implements Chebyshev polynomials with a Chebyshev–Gauss–Lobatto quadrature grid that clusters the grid points near the walls, suitable to resolve the boundary-layer dynamics. Recent studies have shown the capability of the model to solve nonlinear internal wave problems, at laboratory scales, to investigate wave–boundary interaction (Deepwell et al. 2021; Hartharn-Evans et al. 2022), boundary-layer instability (Harnanan et al. 2017; Posada-Bedoya et al. 2024) and sediment resuspension (Olsthoorn & Stastna 2014).

In the present set-up, no-slip and no-flux boundary conditions were imposed on the top and bottom boundaries, with periodic horizontal boundary conditions (figure 2 a). A Chebyshev–Gauss–Lobatto grid was employed in the vertical direction with a clustering of grid points near the top and bottom walls and a uniform grid was used in the horizontal direction. Grid resolution was defined to solve at least 15 points per roughness wavelength in the horizontal direction, and at least 7 Chebyshev grid points below the roughness height, i.e. below $z=h_b$ . These grid resolutions are comparable to those used by Carr et al. (2010) to resolve the corrugated bed, and by Stastna & Lamb (2008) to resolve the flow surrounding bottom topography. The terrain-following coordinates used by the numerical model relax the grid requirements because of the exact definition of the boundary conditions following the variable bottom. Grid resolutions ranged from 2048 $\times$ 256 to 8192 $\times$ 4096 in the horizontal and vertical directions respectively. We conducted additional simulations doubling the grid resolution in both directions for selected cases to verify grid independence at these resolutions (Appendix A). From these, we verified grid independence down to 3 grid points below $z=h_b$ and 8 grid points per roughness wavelength.

2.3. Parameter space

We conducted an extensive parametric study that explored the effect of changes in ISW-induced BBL parameters $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , and roughness parameters $h_b$ and $\lambda _b$ on the nature of the BBL instability. We considered 20 different ISWs over the $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space (figure 1). To investigate the effect of roughness on the BBL instability, we analysed separately the effect of changing the roughness height, while keeping the roughness wavelength constant for each ISW, and vice versa. A total of 96 simulations were carried out for this study: 80 out of the 96 simulations varied the roughness height while keeping the same wavelength for each ISW; the remaining 16 simulations were selected cases to investigate the effect of varying the roughness wavelength.

2.3.1. Roughness height

For each ISW we simulated a flat-bottom case ( $h_b/H=0$ ) and three roughness heights: $h_b/H=10^{-3}$ , $h_b/H=10^{-4}$ and $h_b/H=10^{-5}$ . For a laboratory-scale 1 m deep tank, these correspond to 1, 0.1 and 0.01 mm respectively, which are close to representative roughness scales of typical flume materials (Darby & Chhabra 2017); 0.001 mm–0.01 mm for glass ( $h_b/H=10^{-5}$ ) and 0.025 mm–0.2 mm for smooth concrete ( $h_b/H=10^{-4}$ ). For rough concrete and imperfections, for example, caulking used for joining glass panels, or when considering a sediment bed (e.g. Aghsaee & Boegman 2015; Ghassemi et al. 2022) roughness can be $\sim \mathcal {O}(1)$ mm, represented here by the cases with $h_b/H=10^{-3}$ .

To investigate the effect of roughness height, the roughness wavelength was kept constant for each ISW case. We set the roughness wavelength equal to the wavelength of the most unstable mode of the BBL ( $\lambda _b^{OS}$ ) predicted by linear stability theory (see § 2.4). Therefore, we considered optimal conditions for the roughness elements to favour growth of instabilities in the BBL, as has been shown for surface periodic waves (Blondeaux & Vittori 1994) and SSWs (Scandura 2013). By analysing different ISWs over the $(Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W})$ space, this definition of $\lambda _b$ allowed us to compute the most unstable $(Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W})$ for a given $h_b$ , as often done to define thresholds of BBL instability (e.g.Schlichting 1968; Verschaeve & Pedersen 2014). This was a more consistent approach than, for example, arbitrarily selecting ‘a priori’ a constant roughness wavelength for all ISWs, or considering random roughness, which would make it more difficult to parametrically investigate the effect of roughness wavelength on the BBL instability. Complementary laboratory experiments (unpublished, in preparation) of periodic ISWs propagating over a deformable sediment bed showed that bedforms develop with a characteristic wavelength that matches that of the most unstable mode of the BBL, which also motivated our definition of $\lambda _b$ . In the discussion section, we will show that the conclusions drawn here, based on the parametric definition of the roughness wavelength, remain valid and can be extended to more realistic roughness conditions.

2.3.2. Roughness wavelength

To investigate the effect of roughness wavelength, we repeated the simulations with two additional roughness wavelengths for selected ISWs (cases 3, 8, 13 and 18 in table 1). These additional simulations were only repeated for $h_b/H=10^{-3}$ and $10^{-4}$ . The additional roughness wavelengths were chosen far from the initially simulated wavelength of the most unstable mode: $\lambda _{b}^{-}=\lambda _b^{OS}/3$ and $\lambda _{b}^{+}=3\lambda _b^{OS}$ , considering that $\lambda _{b}^{-}$ is limited by the horizontal grid resolution.

The parameters of the simulated ISWs and their associated BBLs are summarized in table 1. Figure 1 illustrates the selected ISWs in the $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ and $a/H$ versus $Re_w$ parameter space, along with the corresponding experiments by Carr et al. (2008) and Zahedi et al. (2021). The simulations were designed to cover the range of parameters encompassing the various numerical and experimental instability thresholds. This allowed us to investigate the effects of bottom roughness on the different instability parameterizations.

2.4. Space–time linear stability analysis

We conducted a space–time linear stability analysis to: (a) determine the wavelength of the most unstable mode for each ISW, which was set as the roughness wavelength in our simulations (see § 2.3.2), and (b) compute the total linear spatial amplification in the BBL beneath each ISW, which was analysed in conjunction with the thresholds of instability predicted by our simulations.

For each ISW, we analysed the stability properties of the reverse-flow jet within the separated BBL (i.e. upstream of the separation point) in the frame of reference of the progressive ISW. The velocity profiles were uniformly sampled behind the ISW trough, from the 2-D flat-bottom simulations, at a time prior to the onset of instability – if observed. For each velocity profile $U(z)$ , we solved the Orr–Sommerfeld (O–S) equation (Drazin & Reid 1981)

(2.7) \begin{equation} \frac {i}{Re}\left ( \frac {\mathrm {d}^4 \hat {v}}{\mathrm {d} z^4 } - 2\alpha ^2\frac {\mathrm {d}^2 \hat {v}}{\mathrm {d} z^2 } + \alpha ^4 \hat {v} \right )- \left (\alpha U(z) -\omega \right )\left ( \frac {\mathrm {d}^2 \hat {v}}{\mathrm {d} z^2 } -\alpha ^2 \hat {v} \right ) - \alpha \frac {\mathrm {d}^2 U(z)}{\mathrm {d} z^2 }\hat {v}=0, \end{equation}

with boundary conditions

(2.8) \begin{equation} \hat {v}(0)=\frac {{\rm d}\hat {v}}{{\rm d}z}(0)=0, \quad \hat {v}(z\rightarrow \infty )\rightarrow 0 , \quad \frac {{\rm d}\hat {v}}{{\rm d}z}(z\rightarrow \infty )\rightarrow 0. \end{equation}

We determined the set of complex frequencies $\omega =\omega _r+i\omega _i$ associated with a corresponding given set of complex wavenumbers $\alpha =\alpha _r+ i \alpha _i$ for perturbations with the general form

(2.9) \begin{equation} {v}' = \hat {v}(z)e^{i(\alpha x-\omega t)}. \end{equation}

Under this approach, $\omega _r$ and $\alpha _r$ represent the frequency and wavenumber of the linear instability, and $\omega _i$ and - $\alpha _i$ , respectively, represent the growth (time) and amplification (space) rate; where we make the distinction between growth in time and amplification in space. Equation (2.7) was solved using a Chebyshev collocation method (Orszag 1971). Further details of the Orr–Sommerfeld equation solution and its validation are presented in Appendix B.

Figure 3. Snapshots of the non-dimensional vorticity field ( $\omega _y\delta _s/c$ ) after the ISW passage over (a,e) the flat bottom and the rough wall region with three simulated roughness heights: (b,f) $h_b/H=10^{-5}$ , (c,g) $h_b/H=10^{-4}$ , (d,h) $h_b/H=10^{-3}$ . Results for two selected ISWs with $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ : (a,b,c,d) $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=98$ and (e,f,g,h) $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=240$ . Cases correspond to ISW 12 and 14 in table 1, respectively. For each ISW, the roughness wavelength was the wavelength of the most unstable mode of the BBL (see $\lambda _b^{OS}/L_w$ in table 1). Each panel shows a near-bed zoom of the vorticity field over the rough wall region and a zoom of the bottom topography. Note that the domains here are shifted to the right from the schematic in figure 2.

The Orr–Sommerfeld equation assumes steady, incompressible and parallel base flow $U(z)$ . The stability analysis was conducted in the frame of reference of the ISW, where the inviscid flow imposed on the inner layer was steady (e.g. Verschaeve & Pedersen 2014; Posada-Bedoya et al. 2024). The flow in the inner layer is quasi-steady because the rate of change within the inner layer is sufficiently low or close to the reference frame of the solitary wave. Non-parallel effects can be expected to be weak in shallow laminar separated flows (Diwan & Ramesh 2012), with a small aspect ratio (i.e. height-to-length ratio, $h_b/L_b$ ). For this reason, local parallel stability theory works well for laminar separation bubbles (e.g. Marxen et al. 2003). In our simulations, the aspect ratios of streamlines, in the separation region ( $h_b/L_b\approx 0.02-0.04$ ) (e.g. compare $x$ and $z$ scales in figure 2(c)), are within typical values reported in the literature (see figure 3 a in Diwan & Ramesh (2012)) and so we can invoke the arguments of Diwan & Ramesh (2012) to consider a locally parallel BBL, suitable for local linear stability analysis.

We iteratively solved (2.7) for each flat-bottom simulated ISW velocity profile to obtain 2-D maps of ( $\omega _r$ , $\omega _i$ ) as a function of selected sets of ( $\alpha _r$ , $\alpha _i$ ); similar to Jones et al. (2008). The values of $\omega _i$ along ( $\alpha _r$ , $\alpha _i=0$ ) correspond to the solution of the temporal problem of the O-S equation, which predicts the spectrum of growth rates for a set of real wavenumbers $\alpha _r$ . Likewise, the values of $-\alpha _i$ along ( $\omega _r$ , $\omega _i=0$ ) correspond to the solution of the spatial problem of the O-S equation, which predicts the spectrum of amplification rates ( $-\alpha _i$ ) for a set of real frequencies $\omega _r$ . Figure 14 in Appendix C shows an illustrative example of the 2-D map and the corresponding spectra of growth and amplification rates for a selected velocity profile. This procedure was repeated for each velocity profile sampled behind each ISW trough. Note that this spatial-temporal 2-D mapping is often conducted when looking for the presence of convective or absolute instability through the cusp-map technique (Kupfer et al. 1987; Schmid & Henningson 2001).

2.4.1. Most unstable wavelength

For each ISW, we defined a corresponding wavelength representative of the most unstable spatial mode (i.e. largest $-\alpha _i$ ), $\lambda _b^{OS}$ , which was set as the wavelength of the periodic roughness in our 2-D simulations. Given the weak streamwise variability of the BBL flow, the computed most unstable wavelength slightly varied in the streamwise direction. To define a single wavelength representative of the most unstable mode of each BBL, $\lambda _b^{OS}$ was computed as the weighted average of the most unstable wavelengths over the pocket of local instability, with amplification rate as the weighing factor. For selected cases, we verified that the BBL instability, predicted by the 2-D simulations, was not significantly sensitive to variations in $\lambda _b^{OS}$ within $\pm$ 10%; likely due to the continuous spectrum of unstable modes. Values of $\lambda _b^{OS}$ are summarized in table 1.

2.4.2. Amplification factor

We followed the approach of Verschaeve & Pedersen (2014) to compute the spatial amplification of the linear instability beneath the ISW. The amplitude ratio of the instability ( $A=A(x)$ ) relative to the unknown but small value $A_0$ at the separation point $x_{sep}$ (where the region of instability begins, see figure 2(c)) was computed

(2.10) \begin{equation} \log _{10}\left (\frac {A}{A_0} \right )=-\int _{x_{sep}}^{x}\alpha _i(x){\rm d}x. \end{equation}

Here, $-\alpha _i(x)$ is the amplification rate of the most unstable mode computed for each selected velocity profile (i.e. at each location $x$ ). This approach is known as the $N$ -factor method (Herbert 1997), where the $N$ factor is given by the maximum amplitude ratio in the streamwise direction

(2.11) \begin{equation} N=\max \left (\log _{10}\frac {A}{A_0} \right ). \end{equation}

Results of the $N$ -factor method are presented in § 3.3.

3. Results

3.1. Roughness-induced BBL instability

For brevity, we focus the description on selected illustrative ISW cases. A similar dependence of the BBL stability on the roughness parameters was identified for all ISWs in the parameter space.

3.1.1. Effect of roughness height

To show the effects of roughness height, we describe the production of vortex-shedding for ISW 12 ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=98$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ ) and ISW 14 ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=240$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ ). These were chosen because they were both expected to be stable and remained laminar when the ISW passed over the flat-bottom region in our simulations (figure 3 a,e), in agreement with the 2-D simulation-derived threshold formulated by Aghsaee et al. (2012). However, for the $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=240$ case, the BBL was expected to be unstable according to the laboratory experiments by Carr et al. (2008) and Zahedi et al. (2021) (figure 1).

To visualize the effects of varying the roughness height on the production of near-bed vortex shedding, figure 3 shows the vorticity for these cases. The roughness wavelength was equal to that of the most unstable mode of the reverse-flow jet predicted by the O-S equation from the associated flat-bed case (see 2.3.1 and table 1).

For $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=98$ the amplitude of the near-bed vorticity field perturbations increased with increasing roughness height. However, vortex shedding only occurred for $h_b/H=10^{-3}$ ( $h_b/\delta _v=0.44$ ) (figure 3 d). For $h_b/H=10^{-4}$ , the roughness elements introduced instabilities in the near-bed vorticity field, but these were not sufficient to trigger vortex shedding (figure 3 c). Rather, the roughness-induced instabilities decayed after the passage of the ISW. For $h_b/H=10^{-5}$ , the perturbations were much smaller and the BBL flow was nearly insensitive to the roughness, closely resembling the laminar flat-bottom case (cf. panels a and b in figure 3).

For $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=240$ , the BBL was stable over the flat-bottom, but vortex shedding occurred over the rough wall region for the three simulated roughness heights (figure 3 f–h), even for $h_b/H=10^{-5}$ , which was two orders of magnitude smaller than the viscous sublayer thickness: $h_b/\delta _v=0.017$ . This result clearly shows that small-scale roughness, characteristic of smooth flume materials, is a BBL destabilizing mechanism beneath ISWs.

Vortex shedding occurred earlier for larger roughness heights (cf. panels f–h in figure 3). For $h_b/H=10^{-3}$ , vortices detached completely and moved away from the wall over the entire rough wall region. In comparison, for $h_b/H=10^{-4}$ and $h_b/H=10^{-5}$ , the roughness-induced instability was still transitioning after the same elapsed time, with vortices forming near the front end of the rough wall region; thus showing differences in the vertical position of the vortices along the rough bed. Therefore, increasing the roughness height increases the amplitude of the seeding perturbations, thereby facilitating the transition toward vortex shedding.

The nature of the BBL stability does not only respond to changes in bottom roughness but also depends on the ISW-induced BBL parameters $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ and $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ . For example, ISW 12 over $h_b/H=10^{-4}$ (figure 3 c) and ISW 14 over $h_b/H=10^{-5}$ (figure 3 f) had similar values of $h_b/\delta _v$ ; 0.044 and 0.017, respectively. However, the BBL stability was different between these cases, with only ISW 14 triggering vortex shedding. This difference in vortex shedding between different $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ at similar $h_b/\delta _v$ shows that vortex shedding of the BBL, in response to the bottom roughness, is not solely dependent on $h_b/\delta _v$ .

A similar response to changes in roughness height occurred for all waves in the parameter space, with vortex shedding being triggered for a critical $h_b$ that varied depending on $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ and $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , as we discuss later.

3.1.2. Effect of roughness wavelength

To illustrate the effects of roughness wavelength on BBL instability, we considered ISW 3 and ISW 13 (table 1) propagating over roughness with three different wavelengths and the same $h_b/H=10^{-4}$ (see § 2.3.2). The roughness wavelengths ( $\lambda _b^{OS}$ , $\lambda _b^{-}=\lambda _b^{OS}/3$ and $\lambda _b^{+}=3\lambda _b^{OS}$ ) were defined based on spectra of the growth and amplification rates for velocity profiles of the reverse-flow jet within the separated BBL (figure 4) (see § 2.3.2).

Figure 4. (a,e) Snapshots of the near-bed horizontal velocity from simulations of two selected ISWs propagating over a flat bottom. The waves are visualized in the ISW frame of reference, with $x_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0$ at the ISW trough. (b,f) Selected near-bed velocity profiles of the separated BBL in the flat-bottom simulation for the two selected ISWs. Locations of the profiles are indicated by corresponding coloured vertical dashed lines in panels (a) and (e), respectively. (c,g) Growth rate ( $\omega _i$ ) and (d,h) amplification rate ( $-\alpha _i$ ) spectra of unstable modes for each selected velocity profile. Panels show (a,b,c,d) ISW 3 ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=177$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.044$ ) and (e,f,g,h) ISW 13 ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=170$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ ). Vertical dotted lines in panels (c,d,g,h) indicate the roughness wavelength in the three simulated cases.

The simulations (figure 5) support the conclusion that roughness wavelength is a crucial parameter for the stability of the BBL. For the shortest $\lambda _b^-$ and longest $\lambda _b^+$ , vortex shedding did not occur, except for localized instabilities at the edges of the rough wall region (panels a,c and d,f). We attribute those instabilities to the discontinuity in the second derivative of the topography at the transition flat to rough. As the flow at the edges of the roughness region is different from the flow over the periodic bumps, the latter being the focus of the present study, the instabilities at the edges of the rough wall region are inconsequential for the present analysis. Conversely, BBL instability and uniform vortexshedding, over the rough wall region, occurred when the roughness wavelength matched that of the most unstable mode of each BBL (panels b and e). This result suggests that the BBL beneath ISWs behaves as a resonator, with a tendency to preferably amplify perturbations with a wavenumber closer to that of the most unstable mode of the BBL flow. The amplification and vortex shedding are more energetic for cases with larger growth/amplification rates predicted by linear stability theory (cf. panels c,d versus g,h in figure 4). In the discussion section, we present additional simulations using random roughness to further explore the role of roughness wavelength on instability amplification.

Figure 5. Snapshots of the non-dimensional vorticity field ( $\omega _y\delta _s/c$ ) after the ISW passage over the rough wall region for two selected ISWs: (a,b,c) ISW 3A, 3 and 3B ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=177$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.044$ ), respectively, and (d,e,f) ISW 13A, 13 and 13B ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=170$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ ), respectively. All the cases had the same roughness height $h_b/H=10^{-4}$ . For each ISW, three different roughness wavelengths were simulated: (a,d) $\lambda _b^-=\lambda _b^{OS}/3$ , (b,e) $\lambda _b^{OS}/L_w$ , (c,f) $\lambda _b^+=3\lambda _b^{OS}$ (see table 1). Each panel shows a near-bed zoom of the vorticity field over the rough wall region and a zoom of the bottom topography. Note that the domains here are shifted to the right from the schematic in figure 2.

3.2. Bottom boundary layer instability amplification mechanisms

We have shown that changes in both roughness height and wavelength modify the stability of the ISW BBL and provide a path for vortexshedding, in an otherwise stable flow over a flat bottom. We have identified two different vortex-shedding mechanisms for $h_b/\delta _v\approx 1$ and $h_b/\delta _v\approx 10^{-1}$ . Here, we describe the instability amplification mechanisms, in more detail, over the rough wall region for ISW 13 ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=170$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ ). A similar dynamics was observed for other ISWs, but with differing amplification and vortex-shedding rates depending on the set $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $h_b$ and $\lambda _b$ .

For the case of $h_b/\delta _v\approx 1$ , as the ISW passed over the rough wall region, snapshots of the near-bed vorticity and streamlines show local flow separation at the aft of bumps close to the separation point of the BBL (figure 6). The local separation ejected the vortices out of the BBL. This destabilization and vortex shedding occurred early upon the ISW passage and were quite close to the separation point of the BBL, behind the ISW trough. The production of vortex-shedding was similar to the observations by Carr et al. (2010) over corrugated beds, which is more representative of ISW propagation over bottom topography as opposed to small-scale surface roughness.

Figure 6. Snapshots of the non-dimensional vorticity field ( $\omega _y\delta _s/c$ ) and streamlines during ISW passage over the rough wall region for ISW 13 with $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=170$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ . The roughness height $h_b/H=10^{-3}$ ( $h_b/\delta _v=1.03$ ). The vertical dashed line indicates the location of the ISW trough. The star indicates the location of the separation point.

For a significantly lower roughness height, with $h_b/\delta _v\approx 10^{-1}$ , the BBL flow oscillated up and down with the roughness wavelength, without local separation or vortices at the aft of each bump (see figure 7). Rather, the effect of roughness was to seed perturbations in the BBL as fluctuations in the velocity field. The seed perturbations gradually grew with time, within the BBL, until they were of sufficient amplitude for nonlinear effects to become important. After further amplification, the instability field became steeper until streamlines overturned to form vortices that shed from the BBL. Therefore, the seed perturbations forced by small-scale roughness amplified over longer timescales and led to vortex shedding at later times that were farther behind the separation point, compared with $h_b/\delta _v\approx 1$ .

For the smallest $h_b/\delta _v\approx 10^{-2}$ , the same mechanism of roughness-induced noise and BBL amplification occurred, but the timing of the vortex shedding differed due to the smaller magnitude of the seed perturbations on the BBL (not shown). In the absence of roughness, this ISW case remains stable. Therefore, it is expected that with a small enough roughness height, the total BBL amplification of the seed perturbations will be insufficient to overcome viscous damping and produce vortex shedding and the BBL will remain stable.

All cases simulated with $h_b/\delta _v\gtrapprox 1$ led to vortex shedding with a phenomenology similar to that described for figure 6. However, even though for $h_b/\delta _v\lessapprox 10^{-1}$ perturbations were always seeded within the BBL (e.g. figure 7), vortex shedding did not always occur. As described in § 3.1.1, vortex shedding of the ISW-induced BBL over roughness is not only dependent on $h_b/\delta _v$ , but also on the ISW-induced BBL parameters $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ and $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ . This motivated us to investigate a roughness-induced vortex-shedding criterion, as described in the following sections.

Figure 7. Snapshots of the non-dimensional vorticity field ( $\omega _y\delta _s/c$ ) and streamlines during the ISW passage over the rough wall region for ISW 13 with $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=170$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}=0.1$ . The roughness height $h_b/H=10^{-4}$ ( $h_b/\delta _v=0.1$ ). The vertical dashed line indicates the location of the ISW trough. The star indicates the location of the separation point. Panel (a) contains an inset zoom of the bottom roughness elements.

3.3. Amplification ratios

The growth of BBL instabilities depends on how roughness perturbations are amplified in the wake of the ISW. As described in § 2.4.2, we computed the amplification ratio $\log _{10}(A/A_0)$ behind the wave trough for all ISWs from the linear stability analysis using the N-factor method. We only used velocity profiles over the flat bottom, thereby representing the lowest amplification for each ISW case. We plot the amplification ratio behind the wave trough for all simulated ISWs in figure 8. The maximum asymptotic value in each curve corresponds to the amplification factor $N$ , as defined in (2.11).

Figure 8. Amplitude growth curves behind the ISW trough for all simulated ISWs, in the ISW frame of reference, where $x_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ is the streamwise position in the ISW frame of reference, with abscissa zero at the ISW trough. Continuous blue and red dashed lines indicate stable and unstable BBL cases in the flat-bottom simulations respectively. In-line labels identify the ISWs in table 1.

3.4. Vortex-shedding threshold

To summarize the effects of roughness on BBL stability over all of our simulations, we compared $N$ with the experimental and numerical stability thresholds in $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space (from figure 1). The stability threshold for our simulations, over the $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space, was estimated for each roughness height and the flat-bottom case. In all simulations, the roughness wavelength was that of the most unstable mode for each ISW. Therefore, these thresholds must be interpreted as being associated with the most unstable ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ ) for each $h_b$ , as often done to define thresholds of BBL instability (e.g. Schlichting 1968; Verschaeve & Pedersen 2014). However, below, we show that more realistic random roughness also seeds the most unstable wavelength, validating the instability thresholds estimated here.

Figure 9. Stability diagrams in the $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space for all simulated cases. Coloured contours represent the BBL amplification factor ( $N$ ). (a) Flat bottom ( $h_b/H=0$ ), (b) $h_b/H=10^{-5}$ , (c) $h_b/H=10^{-4}$ and (d) $h_b/H=10^{-3}$ . Pink-shaded regions delimited with black dashed lines indicate the approximated threshold regions associated with ranges of $N=N_c$ .

The discrete distribution of our simulated cases over the parameter space (figure 9) requires the computed thresholds of $N$ to be defined as ranges of $N=N_c$ . Hence, these thresholds were visualized as regions delimited by lines of constant $N=N_c$ that separate stable and unstable ISWs simulated for each $h_b/H$ . Threshold regions shifted to the left, i.e. towards smaller $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , and to smaller values of $N=N_c$ , as the roughness height increased. For $h_b/H=10^{-3}$ , the threshold region was roughly vertical, suggesting that the threshold is mostly determined by a critical $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\approx\,75$ , with a weak dependence on $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ . For decreasing $h_b/H=10^{-4}$ , $h_b/H=10^{-5}$ and $h_b/H=0$ , the $N_c$ lines increased in slope, in the $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space, suggesting an increasing dependence on $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ as $h_b/H$ decreased and $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ increased.

In the flat-bottom case limit, only three cases were unstable in our simulations (ISWs 15, 19 and 20), with a critical $N_{c}\approx 3.3-4.0$ (see figure 9 a), which is of similar order as the critical amplification estimated by Verschaeve & Pedersen (2014) (their figure 27) for the numerical simulations of Aghsaee et al. (2012). The observed instability is of the convective type, as previously described by Posada-Bedoya et al. (2024). The threshold region roughly matches the threshold curve proposed by Aghsaee et al. (2012) using 2-D spectral simulations. Discrepancies are likely due to their arbitrary selection of an equation to fit discrete stable and unstable cases (see their figure 12), which differed from our runs and also exhibited some uncertainty through the transition region. Overall, we consider that our 2-D flat-bottomed spectral simulations match the threshold of instability proposed by Aghsaee et al. (2012).

The threshold region for $h_b/H=10^{-5}$ is the closest to the critical threshold ( $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\sim 200$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\sim 0.05-0.07$ ) associated with the laboratory experiments by Carr et al. (2008) and Zahedi et al. (2021), with $N_c\approx 1.5-1.8$ (see figure 9 b). This shows the critical amplification factor to be approximately two orders of magnitude smaller in the laboratory, compared with the 2-D flat-bottomed spectral simulations (both ours and those by Aghsaee et al. (2012)). This is quite close agreement to the estimates by Verschaeve & Pedersen (2014), which also report a difference of two orders of magnitude between the laboratory experiments by Carr et al. (2008) and the simulations by Aghsaee et al. (2012).

4. Discussion

We have shown that bottom roughness changes the stability of the BBL beneath ISWs of depression. The smallest roughness scale considered here, $h_b/H=10^{-5}$ , in some cases two orders of magnitude smaller than the viscous sublayer thickness, can destabilize the BBL and trigger vortex shedding, in an otherwise numerically stable BBL over a flat bottom. This roughness height is of the same order as in typical flume materials used in the laboratory (Darby & Chhabra 2017): 0.01–0.001 mm for glass and 0.2–0.025 mm for smooth concrete. Therefore, the dependence of the BBL stability on the characteristics of the small-scale surface roughness means that glass and smooth concrete should not be considered “smooth” surfaces in a laboratory setting of an ISW-induced BBL flow (see figure 9).

We identified two mechanisms for vortex shedding determined by the ratio $h_b/\delta _v$ . For $h_b/\delta _v\gtrapprox 1$ , vortices were forced directly by the local flow separation at the aft of each bump, more representative of ISW interaction with bottom topography (e.g. Carr et al. 2010). Conversely, for $h_b/\delta _v\lessapprox 10^{-1}$ , the presence of roughness forces the BBL to oscillate with the roughness wavelength, seeding perturbations in the BBL susceptible to being amplified by the BBL flow.

In both cases, amplification occurred preferably for roughness wavelengths close to those of the most unstable mode of the BBL, as predicted by linear instability theory. This resonator-like nature of the BBL flow beneath ISWs has also been reported for the BBL driven by SSWs (Scandura 2013) and periodic monochromatic waves (Blondeaux & Vittori 1994). Despite the similarities between the SSW and ISW problems, it remains a challenge to relate the stability of both BBL flows due to the additional parameters $\Delta \rho /\rho _0$ and $z_{pyc}/H$ that define the ISWs (e.g.Verschaeve & Pedersen 2014).

4.1. Vortex-shedding criterion: random roughness wavelength

We have arbitrarily chosen the most unstable mode of the BBL as the roughness wavelength in our simulations. Given the resonator-like nature of the BBL, this can be considered the most favourable condition for roughness to destabilize the BBL. The thresholds computed above, therefore, correspond to the optimal and most unstable parameters $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ that would destabilize the BBL for each $h_b$ . To test the generality of these results, we performed simulations using more realistic random distributions of roughness wavelength.

Sixteen additional simulations were conducted to validate the generality of the thresholds in § 3.4. The selected cases were chosen from figure 9 immediately above and below the thresholds for $h_b/H=10^{-4}$ and $10^{-5}$ . For each $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , we chose the pair of cases with the largest stable $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ and the smallest unstable $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ , hence eight simulations per $h_b/H$ (figure 10).

Figure 10. Stability diagrams in the $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space for the sixteen cases simulated with random roughness wavelengths for (a) $h_b/H=10^{-5}$ , (b) $h_b/H=10^{-4}$ . Coloured contours represent the BBL amplification factor ( $N$ ). Pink-shaded regions, delimited with black dashed lines, indicate the approximate threshold regions estimated from the sinusoidal roughness simulations (see figure 9).

The random roughness is more representative of realistic beds measured in the laboratory (e.g. Sumer et al. 2010; Scandura 2013; Ghassemi et al. 2022). The rough bed was modelled as the superposition of sinusoids with a wide, nearly continuous spectrum of wavelengths, each with a definite random amplitude, as commonly done for random-roughness numerical studies (e.g. Stastna & Lamb 2008; Scandura 2013). The shortest modelled wavelength was limited by the grid resolution to 8 grid points. The modelled roughness is shown in figure 11(e), along with the band-pass-filtered signal associated with the most unstable wavelength of the BBL (spectral cutoff (0.8–1.2) $\lambda _b^{OS}$ ),which shows regions where the roughness shape follows the most unstable wavelength of the BBL.

To illustrate the BBL interaction with random roughness, figure 11 shows snapshots of the near-bed instantaneous vertical velocity as a selected ISW propagated over the random roughness region. Localized packets of vertical perturbations, behind the separation point, emanated from regions where the roughness had a dominant variability at the wavelength of the most unstable mode of the BBL (figure 11 a–d). The seeded instability packets amplified over time and slowly moved downstream. At subsequent times, new instability packets formed from new regions where the separated BBL interacted with roughness dominated by the most unstable mode (figure 11 a–d). As a consequence, the resulting envelope of the instability packets roughly followed the modulated amplitude of the band-pass-filtered roughness (cf. $w\mid _{\delta _s}/U_2$ versus band-pass-filtered roughness).

In all the sixteen cases, ISW propagation over random roughness predicted the same corresponding stability regimes (i.e. stable versus unstable, see figure 10) as the simulations with sinusoidal roughness. This indicates that the instability thresholds computed for sinusoidal roughness are valid for more realistic random roughness. This confirms that the BBL can resonate with the roughness as long as the random roughness field contains variance at the wavelength of the most unstable mode.

Figure 11. Snapshots of the instantaneous vertical velocity $w/U_2$ (filled contours) during propagation of ISW 13 (table 1) over the random roughness region at (a) $tc/L_w=\,3.9$ , (b) $tc/L_w=\,4.2$ , (c) $tc/L_w=\,4.5$ and (d) $tc/L_w=\,4.8$ . Each panel shows a transect of $w/U_2$ at $z=\delta _s$ (black horizontal dashed line), with the associated scale on the right axis. The dashed lines indicate the upper and lower envelopes of the signal. (e) Random bottom roughness (grey line) and band-pass-filtered roughness (pink) around the wavelength associated with the most unstable mode of the BBL $\lambda _b^{OS}$ (spectral cutoff (0.8–1.2) $\lambda _b^{OS}$ ). Dashed lines indicate the upper and lower envelopes of the band-pass-filtered signal.

4.2. Roughness-induced perturbations as seed noise

For roughness elements completely within the viscous sublayer, we interpret BBL–roughness interaction as a mechanism for seeding perturbations in the BBL. Whether these perturbations introduced by the roughness will be sufficiently amplified to eventually trigger vortex shedding depends on three parameters: (i) the initial amplitude ( $A_0$ ) and (ii) the initial wavelength of the seed perturbations, and (iii) the amplification of the seed perturbations by the BBL flow ( $N$ ). By considering roughness with the most unstable wavelength, we limited the search of critical instability parameters to $A_0$ and $N$ . This is equivalent to computing the most unstable $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ and $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ for each $h_b$ .

Figure 12. Ranges of (a) critical amplification factor ( $N_c$ ) and (b) background noise amplitude ( $A_0$ ) for the different roughness heights ( $h_b/H$ ) and the flat-bottom simulations. Inset in panel (a) shows dispersion on a semilog x-axis, only including $h_b\gt 0$ . Error bars correspond to the ranges of $N_c$ (and associated $A_0$ ) for each threshold region defined in figure 9. Markers are placed in the middle of each range, and used to fit the curves. Shaded regions indicate the approximate range of $N_c$ and noise levels ( $A_0$ ) in the laboratory experiments, the present solver (approximately the same as in the solver of Aghsaee et al. (2012)), and the finite-volume solver of Ellevold & Grue (2023). The dashed line in panel (b) indicates the convergence tolerance of the Generalized Minimum Residual (GMRES) algorithm in the flat-bottom simulations of the present study ( $10^{-7}$ ).

Lines of constant amplification factor in the $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ versus $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ space are consistent with the simulated separation between stable and unstable ISWs, with the threshold shifting to smaller $N_c$ and $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ as $h_b/H$ increased. By considering the most unstable roughness wavelength, the shifting of the instability thresholds to smaller $N_c$ as $h_b/H$ increased is equivalent to increasing the initial amplitude of the seed perturbations (i.e. $A_0$ ), such that smaller amplification is required to reach vortex shedding. Figure 12(a) shows the ranges of critical amplification associated with each threshold as a function of the roughness height. For finite $h_b$ , the critical amplification $N_c$ follows a logarithmic relation

(4.1) \begin{equation} N_c=-0.6579\log _{10}\left ( \frac {h_b}{H}\right ) -1.6718. \end{equation}

By demonstrating that roughness changes the stability characteristics of the BBL, we further hypothesize that our results can be generalized, by associating each roughness height with an equivalent level of seed noise. For the case of a Tollmien–Schlichting wave, the threshold amplitude to trigger secondary instabilities lies at 1% of the free-stream velocity $U_{\infty }$ (Herbert 1988; Verschaeve & Pedersen 2014). Here, we use the same 1% criterion for primary instability to trigger secondary instabilities. We assume $U_{\infty }\sim U_2$ and conduct the analysis based on the actual dimensional values of the simulated $U_2$ . For our ISWs, $U_2\approx 0.1{-}0.2\,\rm m\,s^-{^1}$ for all waves and so we define $U_{\infty }\sim 10^{-1}\,\rm m\,s^-{^1}$ . Assuming that the secondary instability is triggered once its amplitude has grown to $A\approx 0.01U_2$ in all cases, according to (2.10), $A_0=0.01U_210^{-N_c}$ . Hence, we can estimate the background noise level $A_0$ associated with each $N_c$ , which in turn is associated with each roughness height (figure 12 a). Figure 12(b) shows the relation between $h_b/H$ and the background noise level $A_0$ obtained through this procedure, which follows the relation

(4.2) \begin{equation} A_0=0.04697 \left (\frac {h_b}{H} \right )^{0.6579}. \end{equation}

The interpretation of roughness-induced perturbations as noise stems from the description of the BBL–roughness interaction in figure 7, where the BBL is forced to oscillate following the bottom roughness elements, introducing perturbations in the velocity field analogous to noise. From this perspective, for sinusoidal geometry, we can estimate the vertical velocity fluctuations that build the noise as

(4.3) \begin{equation} \text {Noise}\sim {w}'\sim \frac {2h_b}{\lambda _b}u_{jet}, \end{equation}

where $u_{jet}$ is the velocity of the reverse-flow jet where the BBL flow is unstable; roughly estimated as the 20%–30% of $U_2$ from our simulations. We find that ${w}'$ also depends on $\lambda _b$ and $u_{jet}$ . As both parameters vary between ISWs, this would suggest a heterogeneous distribution of seed ${w}'$ and also $A_0$ over the parameter space. However, we verified that $2u_{jet}/\lambda _b\approx 0.3-1.6$ over the parameter space, such that $h_b$ was the dominant factor that determined the order of magnitude of ${w}'$ in (4.3) and, therefore, of the background noise $A_0$ . A detailed analysis of the effect of seed noise on the BBL instability is beyond the scope of this work. However, we have shown that for roughness elements much smaller than the viscous sublayer thickness, the effect of roughness is equivalent to introducing seed noise to the BBL, which is susceptible to being amplified by the BBL flow.

4.3. Critical thresholds for instability: experiments versus 2-D spectral simulations

An immediate consequence of the above analysis is that the definition of general thresholds for instability requires consideration of the background noise level as an additional variable to characterize flow stability. Therefore, to reconcile critical instability thresholds in the laboratory (Carr et al. 2008; Zahedi et al. 2021) and finite-volume 2-D simulations (Thiem et al. 2011; Ellevold & Grue 2023) versus those predicted by 2-D spectral simulations (Diamessis & Redekopp 2006; Aghsaee et al. 2012; Posada-Bedoya et al. 2024), it is necessary to consider the background noise level, in addition to $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ and $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ .

In the laboratory, the lock-release initialization mechanism of the ISW, bottom roughness, sidewall friction, surface waves, instrument placement, etc., are potential sources of additional noise susceptible to being amplified by the BBL. In flat-bottom 2-D simulations, truncation error and initial and boundary conditions are the main noise sources, depending on the accuracy of the numerical solver.

We estimate that unstable BBLs in the laboratory (Carr et al. 2008; Zahedi et al. 2021) have a critical $N_c\approx 1.5{-}1.8$ (figure 9 b), with an associated background noise $A_0\approx 10^{-4.8}{-}10^{-4.5}$ (figure 12). In comparison, in our flat-bottom spectral 2-D simulations, the unstable BBLs have a critical $N_c\approx 3.3{-}4.0$ , with a background noise $A_0\approx 10^{-7}{-}10^{-6.3}$ , which is of the same order as the convergence tolerance of the Generalized Minimum Residual (GMRES) algorithm used for the simulations ( $10^{-7}$ ). Given that our flat-bottom threshold matches that predicted by the spectral 2-D simulations of Aghsaee et al. (2012) (see figure 9 a), the estimated noise is expected to be of the same order in both solvers. Note that Aghsaee et al. (2012) used the spectral solver of Lamb & Nguyen (2009), which was different from ours. This suggests that there is a difference of two orders of magnitude of the background noise between the laboratory experiments versus our 2-D spectral simulations and those by Aghsaee et al. (2012), which explains the discrepancy in instability thresholds. This is consistent with the prediction by Verschaeve & Pedersen (2014), who also estimated a difference of two orders of magnitude in the background noise between the laboratory experiments by Carr et al. (2008) versus the 2-D spectral simulations by Aghsaee et al. (2012). However, their estimated background noise levels assumed $U_{\infty }\sim \mathcal {O}(1)$ , so their actual values for $A$ and $A_0$ differ from our estimates.

Supported by finite-volume 2-D simulations, Ellevold & Grue (2023) proposed a threshold in terms of $a/H$ , $Re_w$ and $d/H$ , which accurately fit the stability threshold from the laboratory experiments by Carr et al. (2008) (see figure 1 b). They attributed the instability to the truncation error of their numerical solver: ‘the growth of the unstable modes arises from the truncation error of the solver at the fifth decimal place’. They defined noise in terms of the vertical velocity fluctuations and quantified it to be $\approx 1{-}2\times 10^{-5}$ (i.e. $\approx 10^{-5}{-}10^{-4.7}$ ). Interestingly, this is very similar to the background noise level $A_0\approx 10^{-4.8}{-}10^{-4.5}$ estimated herein for laboratory experiments, which is within the threshold region associated with the roughness height $h_b/H=10^{-5}$ (figure 9 b). This strongly suggests that the good agreement between the experiments by Carr et al. (2008) and the simulations by Ellevold & Grue (2023) is related to the low accuracy of their finite-volume solver, whose background noise amplitude is comparable to that in the laboratory. In comparison, solvers like SPINS and that used by Aghsaee et al. (2012) have spectral accuracy. According to our previous analysis, and in agreement with the estimates by Verschaeve & Pedersen (2014), the estimated noise level of spectral solvers is approximately two orders of magnitude smaller than in the laboratory (Carr et al. 2008; Zahedi et al. 2021) and the simulations by Ellevold & Grue (2023), explaining the discrepancies between the stability thresholds.

The 1 % criterion for spanwise instabilities invoked above was developed for steady homogeneous boundary-layer flow (Croswell 1985; Herbert 1988). Our 2-D simulations cannot resolve these secondary spanwise instabilities. Even in the quasi-steady reference frame moving with the wave, our boundary layer has spatial variations in the locations of the inflection points, which influence the growth of secondary instabilities (Croswell 1985; Herbert 1988). Despite these limitations, we still used the 1 % criterion for when the flow would transition to turbulence. We note that our conclusions from the analysis above would not change by varying the assumed 1 % criterion for triggering secondary spanwise instability. Changes in magnitude to the assumed percentage value would shift all data points equally along the $A_0$ axis in figure 12(b), while the relative magnitude of the background noise would stay the same. Even a 10 % criterion would not change our interpretation of the results.

As noted by Zahedi et al. (2021), the combined laboratory data suggest the instability threshold is only determined by $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\approx 200$ , independent of $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ . Conversely, our analysis indicates that the instability threshold also depends on $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ . Since the laboratory data only cover a narrow range of $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\approx 0.05{-}0.07$ for $Re_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}\lessapprox 200$ , it is not possible to use those data to unequivocally establish the independence on $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ . Therefore, further experiments are required to validate the predicted dependence of the stability threshold on $P_{{\kern-0.8pt}I{\kern-0.8pt}S{\kern-0.8pt}W}$ .

5. Conclusions

We conducted a parametric study on the effects of roughness on BBL instability beneath ISWs of depression. Our simulations show that small-scale bottom roughness changes the stability of the BBL beneath ISWs, for roughness on the same scale as typical flume materials used in a laboratory setting. This indicates that glass and smooth concrete should not be considered “smooth” surfaces in the laboratory setting of an ISW-induced BBL flow. The BBL beneath ISWs behaves as a resonator that preferentially amplifies perturbations with the wavelength of the most unstable mode, analogous to the BBL beneath periodic surface waves and SSWs interacting with bottom roughness.

We interpret BBL–roughness interaction as a mechanism for forcing seed noise in the BBL. Supported by the $N$ -factor method for Tollmien–Schlichting waves, we propose relations between the roughness height and the seed noise. By considering the background noise level as an additional variable necessary to characterize the flow stability, our results reconcile the discrepancies between critical thresholds of instability in the laboratory and 2-D finite-volume solvers versus those predicted by 2-D spectral simulations.

Our results motivate future research on the interaction of ISW-induced BBL and deformable sediment beds. Specifically, future research should consider how the interaction of the BBL instability driven by bottom roughness can lead to the formation of morphological structures on deformable beds.