2 Reynolds averaged description and turbulent closure

We consider a fluid flow along the $x$ direction over a bed of elevation denoted by $Z$ . Here, $z$ is the crosswise axis normal to the main bed and $y$ is spanwise. Following the standard separation between average quantities and fluctuating ones (denoted by a prime), the equations governing the mean velocity field $u_{i}$ and the pressure $p$ can be written as

(2.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{i}u_{i}=0\quad \text{and}\quad \unicode[STIX]{x1D70C}\unicode[STIX]{x2202}_{t}u_{i}+u_{j}\unicode[STIX]{x2202}_{j}u_{i}=\unicode[STIX]{x2202}_{j}\unicode[STIX]{x1D70F}_{ij}-\unicode[STIX]{x2202}_{i}p,\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{ij}$ contains the Reynolds stress tensor $-\unicode[STIX]{x1D70C}\overline{u_{i}^{\prime }u_{j}^{\prime }}$ . Here, $\unicode[STIX]{x1D70C}$ is the density of the fluid, assumed to be constant. We use a first-order turbulence closure to relate the stress to the velocity gradient. It involves a turbulent viscosity resulting from the product of a mixing length and a mixing frequency, representing the typical eddy length and time scales (Pope Reference Pope2000). The mixing length $\ell$ depends explicitly on the distance from the bed. The mixing frequency is given by the strain rate modulus $|\dot{\unicode[STIX]{x1D6FE}}|=\sqrt{1/2\dot{\unicode[STIX]{x1D6FE}}_{ij}\dot{\unicode[STIX]{x1D6FE}}_{ij}}$ , where we have introduced the strain rate tensor $\dot{\unicode[STIX]{x1D6FE}}_{ij}=\unicode[STIX]{x2202}_{i}u_{j}+\unicode[STIX]{x2202}_{j}u_{i}$ . In a homogeneous situation along the $x$ -axis, the strain rate reduces to $\unicode[STIX]{x2202}_{z}u_{x}$ . For a constant shear stress associated with a shear velocity $u_{\ast }$ , we write

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{xz}=\unicode[STIX]{x1D70C}\ell ^{2}|\unicode[STIX]{x2202}_{z}u_{x}|\unicode[STIX]{x2202}_{z}u_{x}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{z}u_{x}=\unicode[STIX]{x1D70C}|u_{\ast }|u_{\ast },\end{eqnarray}$$

where $\unicode[STIX]{x1D708}$ is the kinematic fluid viscosity. In order to account for both smooth and rough regimes, we adopt here a van Driest-like mixing length (Pope Reference Pope2000),

(2.3) $$\begin{eqnarray}\ell =\unicode[STIX]{x1D705}(z+rd-Z)\left[1-\exp \left(-\frac{(\unicode[STIX]{x1D70F}_{xz}/\unicode[STIX]{x1D70C})^{1/2}(z+sd-Z)}{\unicode[STIX]{x1D708}{\mathcal{R}}_{t}}\right)\right].\end{eqnarray}$$

In this expression, $\unicode[STIX]{x1D705}=0.4$ is the von Kármán constant, $d$ is the sand equivalent bed roughness size and ${\mathcal{R}}_{t}$ is the van Driest transitional Reynolds number, equal to ${\mathcal{R}}_{t}^{0}\simeq 25$ in the homogeneous case of a flat bed (Pope Reference Pope2000). The exponential term suppresses turbulent mixing within the viscous sublayer, close enough to the bed, $z=Z$ . The dimensionless numbers $r=1/30$ and $s=1/3$ are calibrated with measurements of velocity profiles over varied rough walls (Schultz & Flack Reference Schultz and Flack2009; Flack & Schultz Reference Flack and Schultz2010). Here, $rd$ corresponds to the standard Prandtl hydrodynamical roughness extracted by extrapolating the logarithmic law of the wall at vanishing velocity. On the other hand, and this is more original, $sd$ controls the reduction of the viscous layer thickness upon increasing the bed roughness.

Following the work of Hanratty (Reference Hanratty1981), ${\mathcal{R}}_{t}$ cannot be taken as constant but depends on a dimensionless number ${\mathcal{H}}$ that lags behind the pressure gradient following a relaxation equation,

(2.4) $$\begin{eqnarray}a\frac{\unicode[STIX]{x1D708}}{u_{\ast }}\unicode[STIX]{x2202}_{x}{\mathcal{H}}=\frac{\unicode[STIX]{x1D708}}{u_{\ast }^{3}}\unicode[STIX]{x2202}_{x}(\unicode[STIX]{x1D70F}_{xx}-p)-{\mathcal{H}},\end{eqnarray}$$

where $a$ is the multiplicative factor in front of the space lag. We also introduce $b=1/{\mathcal{R}}_{t}^{0}\text{d}{\mathcal{R}}_{t}/\text{d}{\mathcal{H}}>0$ as the relative variation of ${\mathcal{R}}_{t}$ due to the pressure gradient. The values of these empirical parameters have been set to $a=2000$ and $b=35$ , as calibrated on the behaviour of the basal shear stress over a modulated bed (Charru et al. Reference Charru, Andreotti and Claudin2013). Their values control the amplitude and location of the anomaly at the transition from a laminar to a turbulent response with respect to the bed elevation, as illustrated below.

We shall focus here on quasi two-dimensional situations, i.e. on geometries invariant along the $y$ direction. This simplification is a limitation for the study of non-transverse patterns, such as these scallops which have a kind of chevron shape, but is certainly enough to capture the physics of the instability. We can generally write the stress tensor as

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ij}=\unicode[STIX]{x1D70C}(\ell ^{2}|\dot{\unicode[STIX]{x1D6FE}}|+\unicode[STIX]{x1D708})\dot{\unicode[STIX]{x1D6FE}}_{ij}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D712}^{2}\ell ^{2}|\dot{\unicode[STIX]{x1D6FE}}|^{2}\unicode[STIX]{x1D6FF}_{ij},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker symbol (see Fourrière, Claudin & Andreotti (Reference Fourrière, Claudin and Andreotti2010) and references therein). The typical value of the phenomenological constant $\unicode[STIX]{x1D712}$ is in the range $2$ – $3$ , but is of no importance here as we shall see that only the normal stress difference $\unicode[STIX]{x1D70F}_{xx}-\unicode[STIX]{x1D70F}_{zz}$ matters. This closed model allows us to compute the velocity and stress profiles for given boundary conditions, as illustrated in § 4.

Figure 2.

Basal shear stress coefficients (a) ${\mathcal{A}}$ and (b) ${\mathcal{B}}$ computed by the model as functions of $k\unicode[STIX]{x1D708}/u_{\ast }$ . The different colours code for the value of the parameter ${\mathcal{R}}_{d}$ (see legend). The laminar–turbulent transition (around $k\unicode[STIX]{x1D708}/u_{\ast }\simeq 10^{-3}$ ) gradually disappears when ${\mathcal{R}}_{d}$ increases.

Before addressing the linear stability analysis, it is useful to discuss and show the effect of the bed roughness on the hydrodynamical quantities, here taking the shear stress as a typical example. Anticipating § 4, where a sinusoidally modulated bed $Z(x)=\unicode[STIX]{x1D701}\text{e}^{\text{i}kx}$ of wavenumber $k$ is considered, these hydrodynamical equations can be linearised with respect to the small parameter $k\unicode[STIX]{x1D701}$ . In this case, the shear stress is also modulated and correspondingly takes the generic form $\unicode[STIX]{x1D70F}_{xz}=\unicode[STIX]{x1D70C}u_{\ast }^{2}(1+k\unicode[STIX]{x1D701}\text{e}^{\text{i}kx}S_{t})+\text{c.c.}$ , where $S_{t}(\unicode[STIX]{x1D702})$ is a dimensionless function of the rescaled vertical coordinate $\unicode[STIX]{x1D702}=kz$ . We define the two coefficients ${\mathcal{A}}$ and ${\mathcal{B}}$ by $S_{t}(0)={\mathcal{A}}+\text{i}{\mathcal{B}}$ ; these are computed as outputs of the hydrodynamic model. The basal shear stress is then a sinusoidal function of $x$ and these coefficients encode the in-phase and in-quadrature components in response to the bed modulation. The coefficients ${\mathcal{A}}$ and ${\mathcal{B}}$ are functions of $k$ , as displayed in figure 2. The main point, which is the hydrodynamical novelty of this study, is that they are also found to strongly depend on the bed roughness-based Reynolds number ${\mathcal{R}}_{d}=\text{d}u_{\ast }/\unicode[STIX]{x1D708}$ . As described and discussed by Charru et al. (Reference Charru, Andreotti and Claudin2013), these coefficients present, in the smooth limit ( ${\mathcal{R}}_{d}<1$ ), a marked transition between the turbulent regime associated with small wavenumbers $k\unicode[STIX]{x1D708}/u_{\ast }\lesssim 10^{-4}$ and the laminar regime, typically for $k\unicode[STIX]{x1D708}/u_{\ast }\gtrsim 10^{-2}$ . This transition, evidenced by the experimental data of Hanratty and his co-workers (Zilker, Cook & Hanratty Reference Zilker, Cook and Hanratty1977; Frederick & Hanratty Reference Frederick and Hanratty1988), resembles a ‘crisis’, with sharp variations of ${\mathcal{A}}$ and ${\mathcal{B}}$ , allowing, in particular, for negative values of ${\mathcal{B}}$ . The coefficients $a$ and $b$ introduced above have been adjusted to fit these data (Charru et al. Reference Charru, Andreotti and Claudin2013). Crucially, this transition progressively disappears upon increasing the roughness Reynolds number, recovering the reference rough behaviour above ${\mathcal{R}}_{d}\gtrsim 100$ (Fourrière et al. Reference Fourrière, Claudin and Andreotti2010). Beyond the phenomenology, a detailed physical understanding of this laminar–turbulent transition remains a pending hydrodynamical open problem.

3 Scalar transport

We wish now to describe a passive scalar $\unicode[STIX]{x1D719}$ , e.g. the concentration of a chemical species or the temperature, which is transported by the flow. We model its dynamics by a simple advection–diffusion equation,

(3.1) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}+\unicode[STIX]{x2202}_{x}q_{x}+\unicode[STIX]{x2202}_{z}q_{z}=0,\end{eqnarray}$$

where the flux $\boldsymbol{q}$ is the sum of a convective and a diffusive term, $\boldsymbol{q}=\unicode[STIX]{x1D719}\boldsymbol{u}-D\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ . Here, we take a diffusion coefficient proportional to the turbulent viscosity and write

(3.2) $$\begin{eqnarray}D=\frac{\ell ^{2}|\dot{\unicode[STIX]{x1D6FE}}|}{\unicode[STIX]{x1D6FD}_{t}}+\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{t}$ and $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ are the turbulent and viscous Schmidt numbers (or Prandtl numbers for temperature), here taken as constants. A typical value for liquids as well as gasses is $\unicode[STIX]{x1D6FD}_{t}=0.7$ (Gualtieri et al. Reference Gualtieri, Angeloudis, Bombardelli, Jha and Stoesser2017). For the molecular diffusivity, $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ can be estimated from the Stokes–Einstein relation, $\unicode[STIX]{x1D708}/\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}=k_{B}T/6\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}r_{m}$ , where $k_{B}$ is the Boltzmann constant, $T$ is the temperature and $r_{m}$ is the molecular effective radius. The order of magnitude of molecular diffusion for ions or for dissolved CO $_{2}$ in water is $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}=10^{3}$ ; for the diffusion of particles in an ideal gas, a typical value is $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}=1$ .

In the base state for which the bed is homogeneous in $x$ , we can assume that the solid bed dissolves or melts at a small constant velocity and then perform the computation in the moving frame of reference. Equation (3.1) reduces to $\unicode[STIX]{x2202}_{z}q_{z}=0$ , so that the flux is constant,

(3.3) $$\begin{eqnarray}q_{z}=-D\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}=q_{0}.\end{eqnarray}$$

With the above expression (3.2) for $D$ , we obtain the equation for the vertical profile of the scalar field $\unicode[STIX]{x1D719}$ ,

(3.4) $$\begin{eqnarray}\left(\frac{\ell ^{2}|\dot{\unicode[STIX]{x1D6FE}}|}{\unicode[STIX]{x1D6FD}_{t}}+\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}}\right)\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D719}=-q_{0}.\end{eqnarray}$$

To solve this equation, a condition on the bed must be specified. It depends on the nature of the scalar field, but can be formally written in a unique general way. For thermal problems, we impose the bed temperature, $\unicode[STIX]{x1D719}_{0}=T_{0}$ . For dissolution or sublimation problems, we write a Hertz–Knudsen type of law (Eames, Marr & Sabir Reference Eames, Marr and Sabir1997), with a flux of matter that depends on the concentration at the surface,

(3.5) $$\begin{eqnarray}q_{0}=\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D719}_{sat}-\unicode[STIX]{x1D719}_{0}).\end{eqnarray}$$

Because this condition applies on the bed, where the velocity vanishes, no convective contribution to the flux is involved. In the above expression, $\unicode[STIX]{x1D719}_{sat}$ is the concentration at saturation, and $\unicode[STIX]{x1D6FC}$ is a reaction kinetic constant, here expressed as a velocity scale. For simplicity, we take it constant. In the case of sublimation or evaporation, $\unicode[STIX]{x1D6FC}$ is given by the thermal velocity times a desorption probability (Eames et al. Reference Eames, Marr and Sabir1997), leading to $\unicode[STIX]{x1D6FC}$ in the range $1{-}10~\text{m}~\text{s}^{-1}$ . In dissolution problems, however, we expect $\unicode[STIX]{x1D6FC}$ to take much smaller values, of the order of $10^{-5}~\text{m}~\text{s}^{-1}$ (Falter, Atkinson & Coimbra Reference Falter, Atkinson and Coimbra2005). Equation (3.5) can be re-expressed as $\unicode[STIX]{x1D719}_{0}=\unicode[STIX]{x1D719}_{sat}-q_{0}/\unicode[STIX]{x1D6FC}$ , so that in both cases $\unicode[STIX]{x1D719}_{0}$ is given. Finally, as buoyancy effects are usually negligible in such problems, we consider that the momentum equation is not coupled to the scalar field.

4 Linearised equations

We can now proceed to the linear stability analysis, to compute the dispersion relation.

4.1 Definitions and base state

For small enough amplitudes, we can consider a bottom profile of the form

(4.1) $$\begin{eqnarray}Z(x)=\unicode[STIX]{x1D701}\text{e}^{\text{i}kx}\end{eqnarray}$$

without loss of generality (real parts of expressions are understood). Here, $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}/k$ is the wavelength of the bottom and $\unicode[STIX]{x1D701}$ is the amplitude of the corrugation. The case of an arbitrary relief can be deduced by a simple superposition of Fourier modes. We introduce the dimensionless variable $\unicode[STIX]{x1D702}=kz$ , $\unicode[STIX]{x1D702}_{0}=kd$ and the wavenumber-based Reynolds number ${\mathcal{R}}=u_{\ast }/k\unicode[STIX]{x1D708}$ . Primed quantities in this section denote derivatives with respect to $\unicode[STIX]{x1D702}$ . With these notations and following (2.3), the mixing length in the base state can be expressed in a dimensionless form as

(4.2) $$\begin{eqnarray}k\ell \equiv \unicode[STIX]{x1D6F6}(\unicode[STIX]{x1D702})=\unicode[STIX]{x1D705}(\unicode[STIX]{x1D702}+r\unicode[STIX]{x1D702}_{0})\,(1-\exp (-{\mathcal{R}}(\unicode[STIX]{x1D702}+s\unicode[STIX]{x1D702}_{0})/{\mathcal{R}}_{t}^{0})).\end{eqnarray}$$

We define the function ${\mathcal{U}}(\unicode[STIX]{x1D702})$ giving the wind profile in the base state as $u_{x}\equiv u_{\ast }{\mathcal{U}}$ . From (2.2), we see that it must verify the following equation:

(4.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6F6}^{2}|{\mathcal{U}}^{\prime }|{\mathcal{U}}^{\prime }+{\mathcal{R}}^{-1}{\mathcal{U}}^{\prime }=1,\quad \text{or equivalently}\quad {\mathcal{U}}^{\prime }=\frac{-1+\sqrt{1+4\unicode[STIX]{x1D6F6}^{2}{\mathcal{R}}^{2}}}{2\unicode[STIX]{x1D6F6}^{2}{\mathcal{R}}},\end{eqnarray}$$

which must be solved with the boundary condition ${\mathcal{U}}(0)=0$ corresponding to the no-slip condition of the wind at the solid interface. Similarly, we define the function ${\mathcal{P}}(\unicode[STIX]{x1D702})$ for the passive scalar in the base state as $\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{0}\equiv q_{0}{\mathcal{P}}/u_{\ast }$ . From (3.4), it obeys

(4.4) $$\begin{eqnarray}\left(\frac{\unicode[STIX]{x1D6F6}^{2}|{\mathcal{U}}^{\prime }|}{\unicode[STIX]{x1D6FD}_{t}}+\frac{{\mathcal{R}}^{-1}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}}\,\right){\mathcal{P}}^{\prime }+1=0.\end{eqnarray}$$

The boundary condition is ${\mathcal{P}}(0)=0$ , corresponding to the condition $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{0}$ on the bed.

4.2 Linear expansion

The aim of this subsection is to derive a set of closed equations (4.11)–(4.16), which, once integrated with given boundary conditions, provide as outputs the basal stresses, scalar concentration and flux, in order to compute the dispersion relation of the problem. Although technical, the principle of the calculation is simple and relies on the linearisation of the governing equations with respect to the small parameter $k\unicode[STIX]{x1D701}$ , which is proportional to the aspect ratio of the bed modulation. As in Fourrière et al. (Reference Fourrière, Claudin and Andreotti2010), all quantities are generically written as the sum of their homogeneous term (along $x$ ) plus a linear correction generically written as $ck\unicode[STIX]{x1D701}\text{e}^{\text{i}kx}C(\unicode[STIX]{x1D702})$ . The scaling factor $c$ encodes the physical dimension, whereas $C$ is a non-dimensional mode function which describes the vertical profile of the correction. Scaling velocities with $u_{\ast }$ and stresses with $\unicode[STIX]{x1D70C}u_{\ast }^{2}$ , we define the mode functions $U(\unicode[STIX]{x1D702})$ and $W(\unicode[STIX]{x1D702})$ for $u_{x}$ and $u_{z}$ , and the functions $S_{t}(\unicode[STIX]{x1D702})$ , $S_{n}(\unicode[STIX]{x1D702})$ , $S_{xx}(\unicode[STIX]{x1D702})$ and $S_{zz}(\unicode[STIX]{x1D702})$ for $\unicode[STIX]{x1D70F}_{xz}$ , $p-\unicode[STIX]{x1D70F}_{zz}$ , $\unicode[STIX]{x1D70F}_{xx}$ and $\unicode[STIX]{x1D70F}_{zz}$ respectively. For the dimensionless mixing length $k\ell$ , we write its mode function $L$ . The scalar concentration is scaled by $q_{0}/u_{\ast }$ and its mode function is $\unicode[STIX]{x1D6F7}$ ; its flux is scaled by $q_{0}$ and its mode function is defined as $-F$ . Finally, we call $\unicode[STIX]{x1D6E5}$ the mode function of the coefficient of diffusion. These functions are not all independent, and with these notations, we can in particular express the coefficient of diffusion as

(4.5) $$\begin{eqnarray}D=\frac{u_{\ast }}{k}\left[\left(\frac{{\mathcal{U}}^{\prime }\unicode[STIX]{x1D6F6}^{2}}{\unicode[STIX]{x1D6FD}_{t}}+\frac{{\mathcal{R}}^{-1}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}}\right)+\frac{1}{\unicode[STIX]{x1D6FD}_{t}}[(U^{\prime }+\text{i}W)\unicode[STIX]{x1D6F6}^{2}+2\unicode[STIX]{x1D6F6}L{\mathcal{U}}^{\prime }]k\unicode[STIX]{x1D701}\text{e}^{\text{i}kx}\right].\end{eqnarray}$$

Following (2.5), we express the stress functions as follows:

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle S_{t}=({\mathcal{R}}^{-1}+2\unicode[STIX]{x1D6F6}^{2}{\mathcal{U}}^{\prime })(U^{\prime }+\text{i}W)+2\unicode[STIX]{x1D6F6}{{\mathcal{U}}^{\prime }}^{2}L, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle S_{xx}-S_{zz}=4({\mathcal{R}}^{-1}+\unicode[STIX]{x1D6F6}^{2}{\mathcal{U}}^{\prime })\text{i}U=\frac{4\text{i}U}{{\mathcal{U}}^{\prime }}, & \displaystyle\end{eqnarray}$$

where we have used $({\mathcal{R}}^{-1}+\unicode[STIX]{x1D6F6}^{2}{\mathcal{U}}^{\prime })=1/{\mathcal{U}}^{\prime }$ at the zeroth order (equation (4.3)). Linearising the Hanratty equation (2.4), one obtains

(4.8) $$\begin{eqnarray}({\mathcal{R}}+\text{i}a){\mathcal{H}}=\text{i}(S_{xx}-S_{zz}-S_{n})=-\frac{4U}{{\mathcal{U}}^{\prime }}-\text{i}S_{n}.\end{eqnarray}$$

The linear expansion of the scalar equation gives

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle F^{\prime }=\text{i}{\mathcal{P}}U+\left(\text{i}{\mathcal{U}}+\frac{\unicode[STIX]{x1D6F6}^{2}|{\mathcal{U}}^{\prime }|}{\unicode[STIX]{x1D6FD}_{t}}+\frac{{\mathcal{R}}^{-1}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}}\,\right)\unicode[STIX]{x1D6F7}, & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle \text{with}\quad F=\left(\frac{\unicode[STIX]{x1D6F6}^{2}|{\mathcal{U}}^{\prime }|}{\unicode[STIX]{x1D6FD}_{t}}+\frac{{\mathcal{R}}^{-1}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}}\,\right)\unicode[STIX]{x1D6F7}^{\prime }+\frac{1}{\unicode[STIX]{x1D6FD}_{t}}{\mathcal{P}}^{\prime }\left[\unicode[STIX]{x1D6F6}^{2}(U^{\prime }+\text{i}W)+2\unicode[STIX]{x1D6F6}L{\mathcal{U}}^{\prime }\right]-W{\mathcal{P}}.\quad & \displaystyle\end{eqnarray}$$

Combining these equations with the linearised Navier–Stokes equations, we finally obtain six closed equations,

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle U^{\prime }=-\text{i}W+\frac{S_{t}-2\unicode[STIX]{x1D6F6}{{\mathcal{U}}^{\prime }}^{2}L}{{\mathcal{R}}^{-1}+2\unicode[STIX]{x1D6F6}^{2}{\mathcal{U}}^{\prime }}, & \displaystyle\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle W^{\prime }=-\text{i}U, & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle S_{t}^{\prime }=\left(\text{i}{\mathcal{U}}+\frac{4}{{\mathcal{U}}^{\prime }}\right)U+{\mathcal{U}}^{\prime }W+\text{i}S_{n}, & \displaystyle\end{eqnarray}$$

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle S_{n}^{\prime }=-\text{i}{\mathcal{U}}W+\text{i}S_{t}, & \displaystyle\end{eqnarray}$$

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}^{\prime }=\left.\left[F+W{\mathcal{P}}-\frac{{\mathcal{P}}^{\prime }\left(\unicode[STIX]{x1D6F6}^{2}S_{t}+2\unicode[STIX]{x1D6F6}L{\mathcal{U}}^{\prime }\left({\mathcal{R}}^{-1}+\unicode[STIX]{x1D6F6}^{2}{\mathcal{U}}^{\prime }\right)\right)}{\unicode[STIX]{x1D6FD}_{t}\left({\mathcal{R}}^{-1}+2\unicode[STIX]{x1D6F6}^{2}{\mathcal{U}}^{\prime }\right)}\right]\right/\left[\frac{\unicode[STIX]{x1D6F6}^{2}|{\mathcal{U}}^{\prime }|}{\unicode[STIX]{x1D6FD}_{t}}+\frac{{\mathcal{R}}^{-1}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}}\right],\quad & \displaystyle\end{eqnarray}$$

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle F^{\prime }=\text{i}{\mathcal{P}}U+\left(\text{i}{\mathcal{U}}+\frac{\unicode[STIX]{x1D6F6}^{2}|{\mathcal{U}}^{\prime }|}{\unicode[STIX]{x1D6FD}_{t}}+\frac{{\mathcal{R}}^{-1}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}}\right)\unicode[STIX]{x1D6F7}, & \displaystyle\end{eqnarray}$$

where the disturbance to the rescaled mixing length reads

(4.17) $$\begin{eqnarray}\displaystyle L & = & \displaystyle \unicode[STIX]{x1D705}\left\{-1+\text{e}^{-({\mathcal{R}}/{\mathcal{R}}_{t}^{0})(\unicode[STIX]{x1D702}+s\unicode[STIX]{x1D702}_{0})}\left[1-\frac{{\mathcal{R}}}{{\mathcal{R}}_{t}^{0}}(\unicode[STIX]{x1D702}+r\unicode[STIX]{x1D702}_{0})\right.\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left.\frac{{\mathcal{R}}}{{\mathcal{R}}_{t}^{0}}(\unicode[STIX]{x1D702}+s\unicode[STIX]{x1D702}_{0})(\unicode[STIX]{x1D702}+r\unicode[STIX]{x1D702}_{0})\left(\frac{1}{2}S_{t}+\frac{b}{{\mathcal{R}}+\text{i}a}\left(\frac{4U}{{\mathcal{U}}^{\prime }}+\text{i}S_{n}\right)\right)\right]\right\}.\end{eqnarray}$$

4.3 Boundary conditions

Six boundary conditions must be specified to solve the system (4.11)–(4.16), labelled (i)–(vi) below. The upper boundary corresponds to the limit $\unicode[STIX]{x1D702}\rightarrow \infty$ , in which the vertical fluxes of mass and momentum vanish asymptotically. This means that the first-order corrections to the shear stress and to the vertical velocity must tend to zero: (i) $W(\infty )=0$ and (ii) $S_{t}(\infty )=0$ . In practice, we introduce a finite height $H$ (or $\unicode[STIX]{x1D702}_{H}\equiv kH$ ) at which we impose a null vertical velocity and a constant tangential stress $-\unicode[STIX]{x1D70C}u_{\ast }^{2}$ , so that $W(\unicode[STIX]{x1D702}_{H})=0$ and $S_{t}(\unicode[STIX]{x1D702}_{H})=0$ . Then, we consider the limit $H\rightarrow +\infty$ , i.e. when the results become independent of $H$ . On the bed, $z=Z$ , both velocity components must vanish, which gives (iii) $U(0)=-{\mathcal{U}}^{\prime }(0)$ and (iv) $W(0)=0$ . Using (4.3) and (4.2), we can express

(4.18) $$\begin{eqnarray}U(0)=\frac{1-\sqrt{1+4\unicode[STIX]{x1D6F6}(0)^{2}{\mathcal{R}}^{2}}}{2\unicode[STIX]{x1D6F6}(0)^{2}{\mathcal{R}}}=\frac{1-\sqrt{1+4[\unicode[STIX]{x1D705}r\unicode[STIX]{x1D702}_{0}(1-\exp (-s\unicode[STIX]{x1D702}_{0}{\mathcal{R}}/{\mathcal{R}}_{t}^{0}))]^{2}{\mathcal{R}}^{2}}}{2[\unicode[STIX]{x1D705}r\unicode[STIX]{x1D702}_{0}(1-\exp (-s\unicode[STIX]{x1D702}_{0}{\mathcal{R}}/{\mathcal{R}}_{t}^{0}))]^{2}{\mathcal{R}}}.\end{eqnarray}$$

Regarding the passive scalar, we hypothesise that its flux through the upper boundary remains constant, which gives (v) $F(\unicode[STIX]{x1D702}_{H})=0$ . Finally, on the bed, the dissolution-like condition $q_{z}=\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D719}_{sat}-\unicode[STIX]{x1D719})$ , associated with its zeroth order (3.5), leads to (vi) $F(0)=\unicode[STIX]{x1D6FC}/u_{\ast }[{\mathcal{P}}^{\prime }(0)+\unicode[STIX]{x1D6F7}(0)]$ , where ${\mathcal{P}}^{\prime }(0)$ is known from (4.4). Other quantities like $S_{t}(0)$ and $S_{n}(0)$ come as results of the calculation.

As both $F(0)$ and $\unicode[STIX]{x1D719}(0)$ must remain bounded to be compatible with the bulk equations, we expect two simple asymptotic regimes: in the limit of small $\unicode[STIX]{x1D6FC}/u_{\ast }$ , the scalar modulation $\unicode[STIX]{x1D6F7}(0)$ results from the mixing of the base profile and is independent of $\unicode[STIX]{x1D6FC}$ ; the flux modulation $F(0)$ follows and must vanish linearly in $\unicode[STIX]{x1D6FC}/u_{\ast }$ . Conversely, at large $\unicode[STIX]{x1D6FC}/u_{\ast }$ , the situation is opposite: the substrate is so erodible that any disturbance in concentration at the surface would lead to a diverging flux resorbing it. As $\unicode[STIX]{x1D6F7}(0)$ vanishes as $u_{\ast }/\unicode[STIX]{x1D6FC}$ , the flux $F(0)$ results from the mixing of the base state and is therefore constant.

4.4 Interface growth rate

To compute the temporal evolution of the bed elevation, we proceed with $Z(x,t)=\unicode[STIX]{x1D701}\text{e}^{\unicode[STIX]{x1D70E}t+\text{i}\unicode[STIX]{x1D714}t+\text{i}kx}$ , where $\unicode[STIX]{x1D70E}(k)$ is the growth rate and $\unicode[STIX]{x1D714}(k)$ is the angular frequency of the bed pattern along $x$ . The phase propagation speed is therefore $-\unicode[STIX]{x1D714}/k$ with these notations. As the equations are linear in $\unicode[STIX]{x1D719}$ , one can always define the relevant scalar with the appropriate factor so that the evolution equation for the bottom reads $\unicode[STIX]{x2202}_{t}Z=q_{0}-q_{z}(Z)$ . At the linear order, this gives the dispersion relation,

(4.19) $$\begin{eqnarray}\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714}=q_{0}kF(0).\end{eqnarray}$$

5 Results and discussion

The dispersion relation (4.19) is displayed in figure 3, in the limit of small $\unicode[STIX]{x1D6FC}/u_{\ast }$ and at vanishing ${\mathcal{R}}_{d}$ (smooth case). Following the asymptotic behaviour of $F(0)$ in this limit, the relevant rescaling factor for $\unicode[STIX]{x1D70E}$ is $q_{0}\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D708}$ . We see in (a) a range of unstable wavenumbers with a positive growth rate, in which $\unicode[STIX]{x1D70E}$ reaches a maximum value $\unicode[STIX]{x1D70E}_{m}$ at $k_{m}$ . In this range, the propagation velocity changes sign (b), showing that the instability is absolute and not convective. As shown in figure 4(a), the key result is that the unstable band disappears above a critical value ${\mathcal{R}}_{c}$ of the roughness Reynolds number.

Figure 3.

The dispersion relation for the dissolution instability in the limit of small $\unicode[STIX]{x1D6FC}/u_{\ast }$ , for $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}=10^{3}$ and ${\mathcal{R}}_{d}=10$ . (a) Growth rate $\unicode[STIX]{x1D70E}$ . (b) Angular frequency $\unicode[STIX]{x1D714}$ . Solid lines represent positive values. For negative values, $-\unicode[STIX]{x1D70E}$ or $-\unicode[STIX]{x1D714}$ is plotted (dotted lines). The grey region corresponds to the unstable ( $\unicode[STIX]{x1D70E}>0$ ) range of wavenumbers.

Figure 4.

(a) Stability diagram in the plane $({\mathcal{R}}_{d},k\unicode[STIX]{x1D708}/u_{\ast })$ , computed in the limit of small $\unicode[STIX]{x1D6FC}/u_{\ast }$ and for $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}=10^{3}$ . Solid line: marginal stability curve ( $\unicode[STIX]{x1D70E}=0$ ). Grey: unstable zone ( $\unicode[STIX]{x1D70E}>0$ ). Dotted line: location of the most unstable modes ( $\unicode[STIX]{x1D70E}_{m}$ ). (b) Experimental data of Ashton & Kennedy (Reference Ashton and Kennedy1972) on ice ripples: mean wavelength versus inverse of mean flow velocity. The linear fit $\unicode[STIX]{x1D706}\simeq 0.1/U$ (solid line) corresponds to the red arrow in (a). (c) Histogram of the ripple aspect ratio in these experiments. The black arrow corresponds to the aspect ratio giving ${\mathcal{R}}_{d}={\mathcal{R}}_{c}$ .

As illustrated in figure 5, one can understand the instability mechanism as follows. The erosion of the bed is driven by the mass flux $q_{z}$ , itself controlled by the concentration gradient and the coefficient of diffusion (equation (3.3)). The concentration profile, enforced by the base state, is non-homogeneous, decreasing away from the surface. The crests of a modulated bed profile come closer to regions of lower concentration, enhancing the gradient with respect to the surface where $\unicode[STIX]{x1D719}$ is imposed. For a constant $D$ , this peak effect increases the flux and thus the erosion at the crests, and this stabilising situation is what happens at large $k\unicode[STIX]{x1D708}/u_{\ast }$ , when the wavelength is much smaller than the viscous sublayer. When turbulence is dominant, $D$ is not constant any more, but is controlled by turbulent mixing. At small $k\unicode[STIX]{x1D708}/u_{\ast }$ , turbulence is enhanced slightly upstream of the crests, and hence there is stabilising erosion again. For wavenumbers in the intermediate range corresponding to the laminar–turbulent transition, however, turbulence is shifted downstream by means of the adverse pressure gradient (2.4), enhancing mixing and thus erosion in the troughs, which a is destabilising (amplifying) situation. The opposite behaviour of the in-phase modulations of $D$ and $\unicode[STIX]{x1D719}$ is displayed in figure 6 for the whole range of wavenumbers, showing a change of sign corresponding to enhanced mixing in troughs in the presence of the laminar–turbulent transition only.

Figure 5.

(a) Isocontours of the diffusion coefficient above a modulated bed (in white), with a wavenumber close to the most unstable mode, $k\unicode[STIX]{x1D708}/u_{\ast }=5\times 10^{-4}$ . The other parameter values are ${\mathcal{R}}_{d}=0.1$ , $\unicode[STIX]{x1D6FC}/u_{\ast }=0.01$ , $\unicode[STIX]{x1D702}_{H}=10$ and $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}=10^{3}$ . Red (blue) regions correspond to a strong (weak) mixing. (b) Corresponding isocontours of the concentration $\unicode[STIX]{x1D719}$ . Strong mixing leads to a large vertical flux which reduces $\unicode[STIX]{x1D719}$ . (c) The same as (a), but for a smaller wavenumber, $k\unicode[STIX]{x1D708}/u_{\ast }=5\times 10^{-5}$ . (d) Isocontours of the concentration $\unicode[STIX]{x1D719}$ at a larger wavenumber, $k\unicode[STIX]{x1D708}/u_{\ast }=5\times 10^{-2}$ . In all panels, the flow is from left to right.

Figure 6.

In-phase components of the modulation of the diffusion coefficient $\unicode[STIX]{x0394}(0)$ in red and of the basal concentration $\unicode[STIX]{x1D6F7}(0)$ in blue, for ${\mathcal{R}}_{d}=10$ (a) and ${\mathcal{R}}_{d}=100$ (b), as functions of $k\unicode[STIX]{x1D708}/u_{\ast }$ . Dashed lines represent negative values. The grey region is the unstable range of wavenumbers.

Changing the molecular diffusivity, this phenomenology with a range of unstable modes (figure 4 a) remains, but with consequently varying values of ${\mathcal{R}}_{c}$ and $k_{c}$ (figure 7 a). For ${\mathcal{R}}_{c}$ , we clearly identify a low- $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ regime where ${\mathcal{R}}_{c}\propto 1/\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ and a plateau at large $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ . This means that diffusive processes are controlled by whatever is dominant between the diffusion of momentum ( $\unicode[STIX]{x1D708}$ ) and that of dissolved species ( $\unicode[STIX]{x1D708}/\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ ). Moreover, $k_{c}$ is found to be pretty constant, independent of $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ , with a typical value of around $10^{-3}u_{\ast }/\unicode[STIX]{x1D708}$ . This relative invariance is what one can expect for a bandpass instability. Exploring a broad range of the ratio $\unicode[STIX]{x1D6FC}/u_{\ast }$ in figure 7(b), but now rescaling the growth rate by $q_{0}u_{\ast }/\unicode[STIX]{x1D708}$ , we see that $\unicode[STIX]{x1D70E}_{m}$ becomes independent of it, with a crossover at around $\unicode[STIX]{x1D6FC}/u_{\ast }\simeq 10^{-3}$ . These regimes correspond to those discussed for $F(0)$ , in relation to (4.19). Meanwhile, the most unstable wavenumber $k_{m}$ switches from a plateau value to another value, with a small relative change, emphasising again that relevant wavenumbers are fairly insensitive to all parameters.

Figure 7.

(a) Critical value ${\mathcal{R}}_{c}$ , at which the instability disappears, as a function of $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}$ (solid line, left axis), in the limit of small $\unicode[STIX]{x1D6FC}/u_{\ast }$ . Right axis, dotted line: corresponding critical wavenumber $k_{c}$ . (b) Maximum growth rate (solid line, left axis) and corresponding wavenumber (dotted line, right axis) as functions of $\unicode[STIX]{x1D6FC}/u_{\ast }$ , computed for $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D708}}=10^{3}$ and ${\mathcal{R}}_{d}\rightarrow 0$ .

The development of the instability actually increases the bed roughness, suggesting that the pattern eventually selects the wavenumber $k_{c}$ nonlinearly. As a matter of fact, converting the amplitude of the bed perturbation to a sand equivalent roughness size $d$ (Flack & Schultz Reference Flack and Schultz2010), the value ${\mathcal{R}}_{d}\simeq 80$ is reached with a pattern aspect ratio $2\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}$ of the order of 5 % for these values of $k$ , i.e. typically $k\unicode[STIX]{x1D701}\simeq 0.16$ . This value is clearly the upper bound for the linear expansion to make sense, but still reasonably small.

Thomas (Reference Thomas1979) has gathered measurements from various experiments that provide evidence of the global scaling of the selected scallop wavenumber with the viscous length. The fit of these data gives $k\simeq 6\times 10^{-3}u_{\ast }/\unicode[STIX]{x1D708}$ . This is in fair agreement with our results, but to be more quantitative, we must emphasise two difficulties when looking in more detail at specific experiments. A first ambiguity lies in the definition of the wavelength for three-dimensional objects like scallops. For example, in the study of a dissolution pattern on plaster by Blumberg & Curl (Reference Blumberg and Curl1974), $\unicode[STIX]{x1D706}$ is identified as the ratio $\langle \unicode[STIX]{x1D6EC}^{3}\rangle /\langle \unicode[STIX]{x1D6EC}^{2}\rangle$ , where the angle brackets denote an average over measurements of longitudinal scallop sizes $\unicode[STIX]{x1D6EC}$ . With this definition, a selected wavenumber $k\simeq 3\times 10^{-3}u_{\ast }/\unicode[STIX]{x1D708}$ is reported. Another issue in computing the value of $k\unicode[STIX]{x1D708}/u_{\ast }$ from experimental data is the difficulty of relying on an accurate estimate of the shear velocity, as also discussed by Blumberg & Curl (Reference Blumberg and Curl1974). Usually, the mean flow velocity is actually measured, and $u_{\ast }$ is computed assuming a velocity profile, typically the logarithmic law of the wall, which leads to $u_{\ast }\simeq \unicode[STIX]{x1D705}U/\ln H/z_{0}$ , where $H$ is the flow depth and $z_{0}$ is the hydrodynamical roughness, itself related to the bedform amplitude as $z_{0}=rd$ .

Closer to the two-dimensional situation that we consider here, we have more quantitatively investigated the data of Ashton & Kennedy (Reference Ashton and Kennedy1972) for ice ripples. These authors report a clear linear law, $\unicode[STIX]{x1D706}\sim 1/U$ (figure 4 b). They have also measured the ripple aspect ratio. This quantity is widely distributed (figure 4 c), showing a population of emerging bedforms with small aspect ratios, and another one of mature ripples, with an aspect ratio centred around 6 % (figure 4 c). Computing $u_{\ast }$ from $U$ as discussed above, we obtain for these data $k\unicode[STIX]{x1D708}/u_{\ast }\simeq 8\times 10^{-4}$ , in quantitative agreement with the value of $k_{c}$ (figure 4 a). Further experimental studies are needed to investigate the emergence and development of this instability in more detail, and in particular to follow the evolution of the bed roughness over time (Villen, Zheng & Lister Reference Villen, Zheng and Lister2005). Another direction of research is the investigation of fast flows in order to connect scallops generated by water flows with those on meteorites, for which supersonic effects are expected.