2.1 Turbulent flow over a sinusoidal bottom

To investigate the linear stability analysis of the flat bed, we consider a sinusoidal bed perturbation whose elevation profile is of the form

(2.2) $$\begin{eqnarray}Z=\unicode[STIX]{x1D701}\text{e}^{\text{i}k(\cos \unicode[STIX]{x1D6FC}x+\sin \unicode[STIX]{x1D6FC}y-ct)+\unicode[STIX]{x1D70E}t},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ is the angle between the wavevector $\boldsymbol{k}=(k\cos \unicode[STIX]{x1D6FC},k\sin \unicode[STIX]{x1D6FC})$ and the flow direction, as sketched in figure 2(a). The corresponding wavelength is $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}/k$ . The amplitude $\unicode[STIX]{x1D701}$ of the perturbation is such that $k\unicode[STIX]{x1D701}\ll 1$ . Here $\unicode[STIX]{x1D70E}$ is the growth rate and $c$ the propagation velocity of the perturbation. We compute the dispersion relations $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D6FC},k)$ and $c(\unicode[STIX]{x1D6FC},k)$ following Andreotti et al. (Reference Andreotti, Claudin, Devauchelle, Durán and Fourrière2012a ). Denoting by a hat $\hat{\cdot }$ the Fourier transform of the variables, the $x$ - and $y$ -components of the bed shear stress perturbation induced by the undulated bed (2.2) can be written as

(2.3a,b ) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D70F}}_{x}=\unicode[STIX]{x1D70F}({\mathcal{A}}_{x}+\text{i}{\mathcal{B}}_{x})k\hat{Z}\quad \text{and}\quad \hat{\unicode[STIX]{x1D70F}}_{y}=\unicode[STIX]{x1D70F}({\mathcal{A}}_{y}+\text{i}{\mathcal{B}}_{y})k\hat{Z},\end{eqnarray}$$

where $\{{\mathcal{A}}_{x},{\mathcal{A}}_{y}\}$ and $\{{\mathcal{B}}_{x},{\mathcal{B}}_{y}\}$ are the in-phase and in-quadrature components of the stress perturbation with respect to the topography, respectively. These coefficients are functions of the wavenumber, and can be numerically derived from the linearisation of a turbulent model based on the Navier–Stokes equations on a steady sinusoidal bed (see the supplementary material of Andreotti et al. (Reference Andreotti, Claudin, Devauchelle, Durán and Fourrière2012a )). However, this dependence can be neglected for unbounded flows, and it is the case for aeolian dunes which are typically much smaller than the atmospheric boundary layer. For the sake of simplicity, we then use the approximate but analytical expressions:

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{A}}_{x}={\mathcal{A}}_{0}\cos ^{2}\unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}_{x}={\mathcal{B}}_{0}\cos ^{2}\unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{A}}_{y}={\textstyle \frac{1}{2}}{\mathcal{A}}_{0}\cos \unicode[STIX]{x1D6FC}\sin \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}_{y}={\textstyle \frac{1}{2}}{\mathcal{B}}_{0}\cos \unicode[STIX]{x1D6FC}\sin \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

where $\{{\mathcal{A}}_{0},{\mathcal{B}}_{0}\}$ are the coefficients when the dune crests are perpendicular to the flow ( $\unicode[STIX]{x1D6FC}=0^{\circ }$ ). The dependence of these coefficients on $\unicode[STIX]{x1D6FC}$ corresponds to the scaling of the velocity for a potential flow over a sinusoidal topography, i.e. in the outer layer where the fluid inertia is balanced by the pressure gradient (Andreotti et al. Reference Andreotti, Claudin, Devauchelle, Durán and Fourrière2012a ). The factor $1/2$ in the transverse direction is due to the variation of the shear stress as the square of the velocity, which is a common approximation for turbulent flows. Here we use ${\mathcal{A}}_{0}=3.5$ and ${\mathcal{B}}_{0}=2$ , as measured in the field on a flat transverse dune (Claudin et al. Reference Claudin, Wiggs and Andreotti2013).

2.2 Sediment transport

In our model, sediment transport can be characterised by three quantities. The first one is the steady sediment flux on a flat sand bed, also called the saturated flux $q_{sat}$ . It is an increasing function of the basal shear stress $\unicode[STIX]{x1D70F}$ above a threshold $\unicode[STIX]{x1D70F}_{th}$ . As argued by Ungar & Haff (Reference Ungar and Haff1987), the transfer of momentum from the fluid to the saltating grains in the transport layer slows down the flow to the threshold velocity. This results in a linear transport law of the form (Andreotti Reference Andreotti2004; Creyssels et al. Reference Creyssels, Dupont, Ould El Moctar, Valance, Cantat, Jenkins, Pasini and Rasmussen2009; Durán et al. Reference Durán, Claudin and Andreotti2011; Tuan Duc et al. Reference Tuan Duc, Valance, Dupont and Ould El Moctar2011)

(2.8) $$\begin{eqnarray}\boldsymbol{q}_{sat}=\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}_{th})\boldsymbol{t}.\end{eqnarray}$$

Here $\unicode[STIX]{x1D6FA}$ is a dimensional constant, and $\boldsymbol{t}$ is the saturated flux direction (corresponding to the direction of the flow $\boldsymbol{e}_{x}$ for a flat sand bed). In subaqueous environments where bedload is the dominant mode of transport, a power law of exponent $3/2$ is reported for the variation of the sediment flux with the bed shear stress (Meyer-Peter & Müller Reference Meyer-Peter and Müller1948; Bagnold Reference Bagnold1956; Durán, Andreotti & Claudin Reference Durán, Andreotti and Claudin2012). However, the results of the linear stability analysis do not qualitatively depend on the exponent of the transport law.

The second quantity is the transport threshold, which can be expressed in terms of a shear stress $\unicode[STIX]{x1D70F}_{th}$ or a velocity $u_{th}$ defined as $\unicode[STIX]{x1D70F}_{th}\equiv \unicode[STIX]{x1D70C}_{f}u_{th}^{2}$ . For a flat bed with a longitudinal slope $\unicode[STIX]{x2202}_{x}Z$ , this threshold is modified at the linear order in $k\unicode[STIX]{x1D701}$ as $\unicode[STIX]{x1D70F}_{th}(1+\unicode[STIX]{x2202}_{x}Z/\unicode[STIX]{x1D707})$ , where $\unicode[STIX]{x1D707}$ is an effective friction coefficient corresponding to the avalanche slope $\unicode[STIX]{x1D707}\simeq \tan (35^{\circ })=0.7$ . This slope effect has been calibrated for both saltation (Hardisty & Whitehouse Reference Hardisty and Whitehouse1988; Iversen & Rasmussen Reference Iversen and Rasmussen1999) and bedload transport (Fernandez Luque & Van Beek Reference Fernandez Luque and Van Beek1976; Dey Reference Dey2003; Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005).

At the linear order, the transverse slope $\unicode[STIX]{x2202}_{y}Z$ has no impact on the threshold value but affects the transport direction. Following Andreotti, Forterre & Pouliquen (Reference Andreotti, Forterre and Pouliquen2012b ), in analogy with a sliding frictional disc on an inclined surface, particles are driven by the flow but also attracted downwards in the direction of the steepest slope by gravity. This results in a transport direction that reads

(2.9) $$\begin{eqnarray}\boldsymbol{t}=\boldsymbol{e}_{x}+\left(\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D70F}}-\frac{1}{\unicode[STIX]{x1D707}_{y}}\frac{u_{th}}{u_{\ast }}\unicode[STIX]{x2202}_{y}Z\right)\boldsymbol{e}_{y}.\end{eqnarray}$$

This behaviour has been tested for bedload transport by Sekine & Parker (Reference Sekine and Parker1992), showing that the coefficient $\unicode[STIX]{x1D707}_{y}$ is also equal to the avalanche slope $\unicode[STIX]{x1D707}$ . However, saltating grains may have a different dynamics because they are not in permanent contact with the bed. In contrast with the case of longitudinal slope, there are no available experimental data to calibrate (2.9). This expression accounts, however, for the correct symmetry and trend with respect to $\unicode[STIX]{x2202}_{y}Z$ . For the sake of simplicity, we then also take $\unicode[STIX]{x1D707}_{y}=\unicode[STIX]{x1D707}=\tan (35^{\circ })$ , but keep in mind that this may overestimate the impact of the bed transverse slope on the direction of the saltation flux.

The flux deviates from the direction of the base flow $\boldsymbol{e}_{x}$ , due to the $y$ -component of the shear stress induced by the oblique perturbation. At large flow velocity, this deflection is maximum, but it is reduced close to the transport threshold by the slope effect. We thus show lines following the transport direction for different velocity ratios $u_{\ast }/u_{th}$ in figure 2(b). At wind velocities far from the threshold, this slope effect is negligible and the lines follow the shear stress direction. The deviation shortens the upwind path and lengthens the downwind path. Then, the time spent in the downwind side is longer, and the associated deviation prevails, leading to larger deviations on the lee side of the ridges. Close to the transport threshold, the slope affects the transport direction, reducing the drift in the lee zones. The corresponding lines are therefore closer to the base line.

The third and last important feature of sediment transport is related to non-homogeneous or non-steady conditions. The actual flux $\boldsymbol{q}$ at a given location is not determined by the local shear stress (2.3) only because, as soon as there is a change in flow strength or surface property, it needs some space and time to relax towards its new steady value (Sauermann et al. Reference Sauermann, Kroy and Herrmann2001; Andreotti, Claudin & Douady Reference Andreotti, Claudin and Douady2002b ; Durán et al. Reference Durán, Claudin and Andreotti2011; Pähtz et al. Reference Pähtz, Parteli, Kok and Herrmann2014). The response time of the topography is much longer than that of the sediment transport so that temporal equilibrium is assumed between $\boldsymbol{q}$ and $\boldsymbol{q}_{sat}$ . Close to saturation, a first-order linear relaxation process describes the relaxation process, involving a single length scale $L_{sat}$ . In a vectorial form to account for both $x$ - and $y$ -components of the flux, it reads

(2.10) $$\begin{eqnarray}L_{sat}\left(\boldsymbol{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)\boldsymbol{q}=\boldsymbol{q}_{sat}-\boldsymbol{q}.\end{eqnarray}$$

For aeolian saltation, a combination of wind tunnel experiments and field data show that $L_{sat}$ is of the order of a metre and proportional to the grain/fluid density ratio $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}$ , and to the grain size $d$ (Andreotti et al. Reference Andreotti, Claudin and Pouliquen2010).

2.3 Growth rate and most unstable mode

The evolution of the bed elevation profile is governed by the sediment mass conservation equation (also known as the Exner equation in this context):

(2.11) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}Z=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}.\end{eqnarray}$$

Here, the sediment flux $\boldsymbol{q}$ is defined as a volumetric flux, i.e. a volume of grains at the bed packing fraction per unit time and unit length, so that the factor in front of the right-hand side of (2.11) is unity. To deal with dimensionless quantities, we rescale lengths by $L_{sat}$ and times by $L_{sat}^{2}/Q$ , where $Q$ is the reference flux $Q=\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FA}u_{\ast }^{2}$ .

Following Andreotti et al. (Reference Andreotti, Claudin, Devauchelle, Durán and Fourrière2012a ), the dispersion relations can be expressed in a dimensionless form as

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}=\frac{k^{2}}{1+(k\cos \unicode[STIX]{x1D6FC})^{2}}[b_{x}\cos \unicode[STIX]{x1D6FC}+b_{y}\sin \unicode[STIX]{x1D6FC}-k\cos \unicode[STIX]{x1D6FC}(a_{x}\cos \unicode[STIX]{x1D6FC}+a_{y}\sin \unicode[STIX]{x1D6FC})], & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle c=\frac{k}{1+(k\cos \unicode[STIX]{x1D6FC})^{2}}[a_{x}\cos \unicode[STIX]{x1D6FC}+a_{y}\sin \unicode[STIX]{x1D6FC}+k\cos \unicode[STIX]{x1D6FC}(b_{x}\cos \unicode[STIX]{x1D6FC}+b_{y}\sin \unicode[STIX]{x1D6FC})], & \displaystyle\end{eqnarray}$$

where the coefficients are

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle a_{x}={\mathcal{A}}_{x}, & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle b_{x}={\mathcal{B}}_{x}-\cos \unicode[STIX]{x1D6FC}\frac{1}{\unicode[STIX]{x1D707}}\frac{u_{th}^{2}}{u_{\ast }^{2}}, & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle a_{y}=\left(1-\frac{u_{th}^{2}}{u_{\ast }^{2}}\right){\mathcal{A}}_{y}, & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle b_{y}=\left(1-\frac{u_{th}^{2}}{u_{\ast }^{2}}\right)\left({\mathcal{B}}_{y}-\sin \unicode[STIX]{x1D6FC}\frac{1}{\unicode[STIX]{x1D707}}\frac{u_{th}}{u_{\ast }}\right). & \displaystyle\end{eqnarray}$$

Plugging the stress coefficients (2.4)–(2.7) into the expression for the growth rate, we obtain

(2.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}= & & \displaystyle \frac{k^{2}}{1+(\cos \unicode[STIX]{x1D6FC}k)^{2}}\left[\cos \unicode[STIX]{x1D6FC}\left({\mathcal{B}}_{0}-k\cos \unicode[STIX]{x1D6FC}{\mathcal{A}}_{0}\right)\left(\cos ^{2}\unicode[STIX]{x1D6FC}+\frac{r}{2}\sin ^{2}\unicode[STIX]{x1D6FC}\right)\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{1}{\unicode[STIX]{x1D707}}\frac{u_{th}}{u_{\ast }}\left(\frac{u_{th}}{u_{\ast }}\cos ^{2}\unicode[STIX]{x1D6FC}+r\sin ^{2}\unicode[STIX]{x1D6FC}\right)\right],\end{eqnarray}$$

where $r\equiv 1-(u_{th}/u_{\ast })^{2}$ quantifies the transport in the direction transverse to the flow (equations (2.9), (2.16) and (2.17)). This function $\unicode[STIX]{x1D70E}$ is shown in figure 3 in the plane $(\unicode[STIX]{x1D6FC},k)$ for the velocity ratios $u_{\ast }/u_{th}=1.6$ and $u_{\ast }/u_{th}=2.5$ . The growth rate is positive (unstable modes) for small wavenumbers (large wavelengths) and transverse flows ( $\unicode[STIX]{x1D6FC}\rightarrow 0$ ). Large wavenumbers are stabilised by the flux relaxation process; longitudinal modes are stabilised by the slope effect on the transport direction. The marginal mode ( $\unicode[STIX]{x1D70E}=0$ ) corresponds to the cutoff wavenumber

(2.19) $$\begin{eqnarray}k_{c}(\unicode[STIX]{x1D6FC})=\frac{1}{\cos \unicode[STIX]{x1D6FC}\,{\mathcal{A}}_{0}}\left[{\mathcal{B}}_{0}-\frac{\displaystyle \frac{u_{th}}{u_{\ast }}\left(\displaystyle \frac{u_{th}}{u_{\ast }}\cos ^{2}\unicode[STIX]{x1D6FC}+r\sin ^{2}\unicode[STIX]{x1D6FC}\right)}{\unicode[STIX]{x1D707}\left(\cos ^{2}\unicode[STIX]{x1D6FC}+\displaystyle \frac{r}{2}\sin ^{2}\unicode[STIX]{x1D6FC}\right)}\right].\end{eqnarray}$$

The selected pattern is the one that maximises the growth rate, and we denote its wavenumber by $k_{max}$ and its orientation by $\unicode[STIX]{x1D6FC}_{max}$ . Taking the derivative of the growth rate (2.12) with respect to $k$ at constant $\unicode[STIX]{x1D6FC}$ , the most unstable wavenumber $k_{\unicode[STIX]{x1D6FC}}$ for a given dune orientation reads

(2.20) $$\begin{eqnarray}k_{\unicode[STIX]{x1D6FC}}=\frac{1}{\cos \unicode[STIX]{x1D6FC}}(X^{1/3}-X^{-1/3}),\end{eqnarray}$$

where

(2.21) $$\begin{eqnarray}X=\sqrt{1+[\cos \unicode[STIX]{x1D6FC}\,k_{c}(\unicode[STIX]{x1D6FC})]^{2}}+\cos \unicode[STIX]{x1D6FC}\,k_{c}(\unicode[STIX]{x1D6FC}).\end{eqnarray}$$

As the saturation length is significantly smaller than the dune wavelength, we can expand the above expressions in the limit $k_{c}(\unicode[STIX]{x1D6FC})\ll 1$ , and we obtain

(2.22) $$\begin{eqnarray}k_{\unicode[STIX]{x1D6FC}}\sim {\textstyle \frac{2}{3}}k_{c}(\unicode[STIX]{x1D6FC}).\end{eqnarray}$$

The overall unstable mode is then $k_{max}=k_{\unicode[STIX]{x1D6FC}_{max}}$ . Under a unidirectional flow, the selected orientation is perpendicular to the flow, i.e $\unicode[STIX]{x1D6FC}_{max}=0^{\circ }$ , which gives

(2.23) $$\begin{eqnarray}k_{max}\sim \frac{2}{3}k_{c}(0)=\frac{2}{3{\mathcal{A}}_{0}}\left({\mathcal{B}}_{0}-\frac{1}{\unicode[STIX]{x1D707}}\frac{u_{th}^{2}}{u_{\ast }^{2}}\right),\end{eqnarray}$$

in agreement with the results of the two-dimensional linear stability analysis (Andreotti et al. Reference Andreotti, Claudin and Douady2002a ; Fourrière et al. Reference Fourrière, Claudin and Andreotti2010; Charru et al. Reference Charru, Andreotti and Claudin2013).

Figure 3.

Growth rate $\unicode[STIX]{x1D70E}$ as a function of the wavenumber $k$ and the orientation angle $\unicode[STIX]{x1D6FC}$ of the bed perturbation, for $u_{\ast }/u_{th}=1.6$ (a,c) and $u_{\ast }/u_{th}=2.5$ (b,d). (a,b) Show the reference transverse case $\unicode[STIX]{x1D6FC}=0^{\circ }$ . (c,d) Show the growth rate in the plane $(\unicode[STIX]{x1D6FC},k)$ with contour lines as indicated. Dashed line: location of the maximum growth rate for each $\unicode[STIX]{x1D6FC}$ . Black dot: global maximum growth rate, corresponding to the most unstable mode $(\unicode[STIX]{x1D6FC}_{max},k_{max})$ . Here, for a unidirectional wind regime, $\unicode[STIX]{x1D6FC}_{max}=0^{\circ }$ , independently of $u_{\ast }/u_{th}$ .

Importantly, we see that the unstable mode exists only if the velocity ratio $u_{\ast }/u_{th}$ is larger than $1/(\sqrt{{\mathcal{B}}_{0}\unicode[STIX]{x1D707}})$ . This condition is always satisfied with the above values of ${\mathcal{B}}_{0}$ and $\unicode[STIX]{x1D707}$ because $1/(\sqrt{{\mathcal{B}}_{0}\unicode[STIX]{x1D707}})\simeq 0.8<1\leqslant u_{\ast }/u_{th}$ . Another related effect is that the selected wavelength $\unicode[STIX]{x1D706}_{max}=2\unicode[STIX]{x03C0}/k_{max}$ tends to increase close to the threshold $u_{\ast }/u_{th}\rightarrow 1$ . For infinitely large values of this velocity ratio, the slope effect on the transport direction vanishes and suppresses the stabilisation of longitudinal modes, so that $\unicode[STIX]{x1D70E}$ diverges in the double limit $k\rightarrow +\infty$ and $\unicode[STIX]{x1D6FC}\rightarrow \pm 90^{\circ }$ . Numerically, for the adopted values of ${\mathcal{A}}_{0}$ , ${\mathcal{B}}_{0}$ and $\unicode[STIX]{x1D707}$ , the regular shape of the growth rate, with a well-defined global maximum, is lost above $u_{\ast }/u_{th}\gtrsim 15$ (an unrealistic regime in practice).

Figure 4.

Characteristics of a bidirectional flow regime. The sketch shows the geometry of the system. The $x$ -axis is defined as the bisector of the divergence angle $\unicode[STIX]{x1D703}$ between the two wind directions (black arrows, unit vectors $\boldsymbol{e}_{1,2}$ ). For the linear stability analysis, we consider a sinusoidal mode characterised by a wavenumber $k$ and an angle $\unicode[STIX]{x1D6FC}$ with respect to the transverse direction ( $y$ -axis). The grey arrows are the directions of the resultant sand flux $\boldsymbol{Q}_{\unicode[STIX]{x1D6F4}}$ corresponding to different values of the transport ratio $N=Q_{1}t_{1}/Q_{2}t_{2}$ . Inset (a): alternation between the two wind directions with respect to time. Inset (b): variation of the amplitude of the resultant flux $Q_{\unicode[STIX]{x1D6F4}}$ for $Q_{1}=Q_{2}$ as a function of $\unicode[STIX]{x1D703}$ . From dark to light greys: increasing values of $N\in \{1,2,3\}$ .

3.1 Linear contribution to the growth rate

We consider two wind directions making an angle $\unicode[STIX]{x1D703}$ , called the divergence angle, whose bisector defines the $x$ -axis (figure 4). The wind blows alternately in directions of the unitary vectors $\boldsymbol{e}_{1}$ and $\boldsymbol{e}_{2}$ over durations $t_{1}$ and $t_{2}$ , and with reference fluxes $Q_{1}$ and $Q_{2}$ . In this bidirectional flow regime, we define the resultant flux $\boldsymbol{Q}_{\unicode[STIX]{x1D6F4}}$ as

(3.1) $$\begin{eqnarray}\boldsymbol{Q}_{\unicode[STIX]{x1D6F4}}=\frac{t_{1}}{t_{1}+t_{2}}Q_{1}\boldsymbol{e}_{1}+\frac{t_{2}}{t_{1}+t_{2}}Q_{2}\boldsymbol{e}_{2}.\end{eqnarray}$$

The intensity of this flux decreases with $\unicode[STIX]{x1D703}$ as the two winds compensate each other (figure 4 b). We also define a time-averaged reference flux as

(3.2) $$\begin{eqnarray}Q=\frac{t_{1}}{t_{1}+t_{2}}Q_{1}+\frac{t_{2}}{t_{1}+t_{2}}Q_{2}.\end{eqnarray}$$

The period of wind reorientation $t_{1}+t_{2}$ is assumed to be small in comparison to the characteristic growth time of the bed perturbation of the order of $(1/\unicode[STIX]{x1D70E}_{max})L_{sat}^{2}/Q$ (Rubin & Ikeda Reference Rubin and Ikeda1990). Thus, the flow changes direction many times while the dune is still in its linear growth regime, during which its orientation and wavelength remain constant. The growth rate $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6F4}}$ of a mode characterised by $(\unicode[STIX]{x1D6FC},k)$ can then be expressed as the sum of the contributions of each wind:

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6F4}}=\frac{1}{N+1}[N\unicode[STIX]{x1D70E}(Q_{1},\unicode[STIX]{x1D6FC}_{1},k)+\unicode[STIX]{x1D70E}(Q_{2},\unicode[STIX]{x1D6FC}_{2},k)],\end{eqnarray}$$

where $N=Q_{1}t_{1}/Q_{2}t_{2}$ is the transport ratio between the two winds (Rubin & Hunter Reference Rubin and Hunter1987; Rubin & Ikeda Reference Rubin and Ikeda1990; Reffet et al. Reference Reffet, Courrech du Pont, Hersen and Douady2010; Courrech du Pont et al. Reference Courrech du Pont, Narteau and Gao2014). The growth rates $\unicode[STIX]{x1D70E}(Q_{1,2},\unicode[STIX]{x1D6FC}_{1,2},k)$ are computed from that derived under a unidirectional flow (2.18) using

(3.4) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{1,2}=\left[\left(\unicode[STIX]{x1D6FC}\pm \frac{\unicode[STIX]{x1D703}}{2}+90^{\circ }\right)\!\!\!\hspace{0.6em}{\rm mod}\hspace{0.2em}180^{\circ }\right]-90^{\circ },\end{eqnarray}$$

which are the angles between the pattern and the perpendicular to each flow direction (figures 2 and 4). For simplicity, we further assume that the winds blow with the same shear velocity $u_{\ast }$ , so that $Q_{1}=Q_{2}=Q$ . As a consequence, different values of $N$ correspond to different blowing times.

The growth rate $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6F4}}$ and the location of its maximum $(\unicode[STIX]{x1D6FC}_{max},k_{max})$ are shown in figure 5 for some values of $N$ and $\unicode[STIX]{x1D703}$ . This maximum integrates the effect of the two winds and gives the orientation and wavenumber of the emerging dune pattern. These selected dune orientations and wavelengths are the main output of the analysis. A third output is the propagation velocity of the selected mode. In the next subsections, we quantitatively discuss the behaviour of these selected quantities with respect to $N$ , $\unicode[STIX]{x1D703}$ and $u_{\ast }/u_{th}$ .

Figure 5.

Growth rate $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6F4}}$ as a function of the wavenumber $k$ and the orientation $\unicode[STIX]{x1D6FC}$ of the bed perturbation for $u_{\ast }/u_{th}=5$ and different values of the divergence angle $\unicode[STIX]{x1D703}$ (columns) and the transport ratio $N$ (rows). Black dots: most unstable mode $(\unicode[STIX]{x1D6FC}_{max},k_{max})$ .

3.2 Selection of the pattern orientation

As shown in figure 6, the orientation of the most unstable mode changes from transverse ( $\unicode[STIX]{x1D6FC}_{max}\approx 0^{\circ }$ ) when the divergence angle $\unicode[STIX]{x1D703}$ is small, to longitudinal ( $\unicode[STIX]{x1D6FC}_{max}\approx -90^{\circ }$ ) when the winds tends to be opposite (large divergence angle). This transition is the signature of a change in the normal to the crest fluxes (Rubin & Hunter Reference Rubin and Hunter1987; Courrech du Pont et al. Reference Courrech du Pont, Narteau and Gao2014). When the divergence angle is smaller than $90^{\circ }$ , the crests are more likely to be perpendicular to the resulting flux, as the sand is mainly transported along the $x$ -direction. By contrast, when the divergence angle is larger than $90^{\circ }$ , the crests are more likely to be aligned with the resulting flow, as the sand is mostly transported along the $y$ -direction. As discussed below, this transition depends on the transport ratio $N$ and the velocity ratio $u_{\ast }/u_{th}$ .

Let us first consider a fixed wind strength ( $u_{\ast }/u_{th}=1.6$ in figure 6 a) and discuss the effect of $N$ . For winds of the same transport capacity ( $N=1$ ), the switch occurs at the transition angle $\unicode[STIX]{x1D703}=90^{\circ }$ , with a discontinuous jump from a pattern perpendicular to the resultant flux direction ( $\unicode[STIX]{x1D6FC}_{max}=0^{\circ }$ ) to a pattern aligned with the resultant flux direction ( $\unicode[STIX]{x1D6FC}_{max}=-90^{\circ }$ ). This transition from a purely transverse to a purely longitudinal dune orientation is specific to the case $N=1$ . When one of the two winds transports more than the other ( $N>1$ ), intermediate orientations (oblique patterns) more perpendicular to the prevailing wind are permitted. There is still a discontinuous transition at the transition angle, but the jump decreases in amplitude as $N$ increases. This discontinuity is associated with two distinct maxima of equal growth rates (figure 5, central column). In the limit of unidirectional flows ( $N\rightarrow +\infty$ ), the variation of $\unicode[STIX]{x1D6FC}_{max}$ tends to a continuous linear function, which simply corresponds to the orientation perpendicular to the dominant wind.

Figure 6.

Pattern orientation $\unicode[STIX]{x1D6FC}_{max}$ as a function of the divergence angle $\unicode[STIX]{x1D703}$ for various transport $N$ and velocity $u_{\ast }/u_{th}$ ratios. (a) Solid lines: predictions of the linear stability analysis for $u_{\ast }/u_{th}=1.6$ . Dashed lines: predictions of Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014) for $\unicode[STIX]{x1D6FE}=1.7$ . From dark to light greys: increasing values of $N\in \{1,2,3,+\infty \}$ . (b–d) Solid black lines: reference curves shown in (a) for $u_{\ast }/u_{th}=1.6$ . Dashed and dotted black lines: predictions of the linear stability analysis for $u_{\ast }/u_{th}=5$ and $u_{\ast }/u_{th}=1.2$ , respectively. Dashed and dotted white lines: predictions of Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014) for $\unicode[STIX]{x1D6FE}=0.5$ and $\unicode[STIX]{x1D6FE}=5$ , respectively. Light grey region: all possible orientations predicted by this analysis when $u_{\ast }/u_{th}$ is increased from $1.005$ to $15$ . Dark grey region: all possible orientations predicted by Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014) when $\unicode[STIX]{x1D6FE}$ is varied from $0$ to $+\infty$ (3.5). Note that, in the strict case of $N=1$ , the growth rate has an equal maximum for $(-\unicode[STIX]{x1D6FC}_{max},k_{max})$ .

The effect of the velocity ratio $u_{\ast }/u_{th}$ , displayed in figure 6 with light grey regions, can be summarised as follows. Strong winds tend to enhance the orientation discontinuity at $\unicode[STIX]{x1D703}=90^{\circ }$ , keeping $\unicode[STIX]{x1D6FC}_{max}$ either closer to $0^{\circ }$ ( $\unicode[STIX]{x1D703}<90^{\circ }$ ) or to $-90^{\circ }$ ( $\unicode[STIX]{x1D703}>90^{\circ }$ ). This corresponds to a selected orientation less perpendicular to the prevailing wind. By contrast, winds close to the threshold induce a more continuous change in orientation, i.e. closer to the perpendicular to the prevailing wind. When $u_{\ast }/u_{th}\gg 1$ , slope effects become negligible and oblique modes stabilised by transverse transport are less damped for large velocity ratios $u_{\ast }/u_{th}$ (figure 3 in the case of a unidirectional wind). When the growth rates associated with the two winds are combined, the ‘cost’ to remain away from the perpendicular to the prevailing wind is not so strong, and the selected orientation reflects this trend. When $u_{\ast }/u_{th}\rightarrow 1$ , slope effects are strong, and mainly affect longitudinal transport, i.e. the coefficient $b_{x}$ (2.15), where two terms compete. The first one, of hydrodynamic origin, scales with $\cos ^{2}\unicode[STIX]{x1D6FC}_{1,2}$ , as deduced from the properties of the potential flow. The second one is related to the variation of the threshold with the slope, i.e. proportional to $-\text{cos}\,\unicode[STIX]{x1D6FC}_{1,2}$ . This difference in the dependence with $\unicode[STIX]{x1D6FC}_{1,2}$ favours a crest orientation perpendicular to the prevailing wind, for which the destabilising hydrodynamics is the strongest.

3.3 Comparison with the dimensional analysis of Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014)

The prediction of the dune orientation by this linear stability analysis can be directly compared to that of Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014). Based on an estimation of the sand flux at the dune crest, these authors computed by dimensional analysis the dune growth rate accounting for the wind speed-up by the topography (appendix A). The selected pattern orientation is then such that it maximises

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6F4}}\propto N\cos \unicode[STIX]{x1D6FC}_{1}+N\unicode[STIX]{x1D6FE}\cos ^{2}\unicode[STIX]{x1D6FC}_{1}+\cos \unicode[STIX]{x1D6FC}_{2}+\unicode[STIX]{x1D6FE}\cos ^{2}\unicode[STIX]{x1D6FC}_{2},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ quantifies the relative increase of sand flux between the flat bed and the dune crest. For $\unicode[STIX]{x1D6FE}=1.7$ (a value also derived from field observations and numerical models by Gao et al. (Reference Gao, Narteau, Rozier and Courrech Du Pont2015b )), the results of that model are found to quantitatively match our predictions for $u_{\ast }/u_{th}=1.6$ (see figure 6). When $\unicode[STIX]{x1D6FE}$ is varied, the curve $\unicode[STIX]{x1D6FC}_{max}(\unicode[STIX]{x1D703})$ changes, and the dark grey area in figure 6(b–d) displays the whole range of possible orientations, showing an overall similar behaviour.

The present linear analysis corresponds to the limit $k\unicode[STIX]{x1D701}\ll 1$ , whereas the computation of Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014) considers a dune with a finite aspect ratio, for which nonlinearities are expected. In a similar fashion to Courrech du Pont (Reference Courrech du Pont2015) and Gao et al. (Reference Gao, Narteau, Rozier and Courrech Du Pont2015b ) (in supplementary material), the parameter $\unicode[STIX]{x1D6FE}$ can be related to the parameters of the present model as

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\frac{|\hat{q}_{sat}|}{q_{sat}}=\frac{k\unicode[STIX]{x1D701}}{r}\sqrt{{\mathcal{A}}_{0}^{2}+{\mathcal{B}}_{0}^{2}},\end{eqnarray}$$

where $\hat{q}_{sat}$ is the Fourier transform of the perturbation of the saturated flux. With the values ${\mathcal{A}}_{0}=3.5$ , ${\mathcal{B}}_{0}=2$ and $u_{\ast }/u_{th}=1.6$ (or equivalently $r\simeq 0.7$ , see (2.18)), the parameter $\unicode[STIX]{x1D6FE}=1.7$ corresponds to a bed perturbation of aspect ratio of $2\unicode[STIX]{x1D701}/\unicode[STIX]{x1D706}=1/12$ . This value is consistent with the typical aspect ratio of fully developed dunes (Parteli et al. Reference Parteli, Schwämmle, Herrmann, Monteiro and Maia2006; Baddock, Livingstone & Wiggs Reference Baddock, Livingstone and Wiggs2007; Elbelrhiti, Andreotti & Claudin Reference Elbelrhiti, Andreotti and Claudin2008). This general consistency between the two approaches emphasises the role played by the normal-to-the-crest components of transport in the selection of dune orientation during the linear but also nonlinear phase of dune growth.

3.4 Selection of the pattern wavenumber

The second key output of the linear stability analysis is the selected pattern wavenumber $k_{max}$ . Remember that $k$ is here made dimensionless by $L_{sat}$ , the only characteristic length scale of the problem. The shape of $k_{max}$ as a function of the divergence angle $\unicode[STIX]{x1D703}$ is displayed in figure 7, for different values of the transport $N$ and velocity $u_{\ast }/u_{th}$ ratios. We can see a non-monotonic behaviour, with variations that can reach $\simeq 30\,\%$ close to the transition angle $\unicode[STIX]{x1D703}=90^{\circ }$ . Interestingly, these variations can be positive or negative, depending on whether the wind strength is far from or close to the transport threshold.

For symmetric winds ( $N=1$ ) and in the limit of large shear velocities, a geometric argument sketched in the insets of figure 7(a) can explain the dependence of the wavelength on the divergence angle (solid curves in that figure). Let us call $\unicode[STIX]{x1D706}_{0}$ the selected wavelength in the unidirectional regime, and $k_{0}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{0}$ the corresponding wavenumber. Far from the transport threshold, and thus without any slope effect, the balance between stabilising and destabilising processes occurs along the wind direction. In a bidirectional wind regime, the selected orientation is oblique to each flow direction, and we can infer that the growth rate is maximum when the wavelength $\unicode[STIX]{x1D706}_{0}$ is selected in the direction of each wind.

Figure 7.

Selected wavenumber $k_{max}$ as a function of the divergence angle $\unicode[STIX]{x1D703}$ for various transport $N$ and velocity $u_{\ast }/u_{th}$ ratios. For comparison with figure 11, note that for $\unicode[STIX]{x1D706}_{max}=2\unicode[STIX]{x03C0}/k_{max}$ the increasing/decreasing trends are opposite. (a) Solid lines: prediction of the linear stability analysis for $u_{\ast }/u_{th}=5$ . Dashed lines: prediction using the geometric scalings, sketched in the insets for $N=1$ (equations (3.12) and (3.13)). From dark to light greys: increasing values of $N\in \{1,2,3,+\infty \}$ . (b) Prediction of the linear stability analysis for $N=1$ . From dark to light greys: increasing values of $u_{\ast }/u_{th}\in \{1.6,2.5,5\}$ .

This argument can be made more precise as follows. Far from the transport threshold, the second term in (2.19) associated with the slope effects vanishes asymptotically as $u_{th}/u_{\ast }\rightarrow 0$ , and this leads to

(3.7) $$\begin{eqnarray}k_{\unicode[STIX]{x1D6FC}}\sim \frac{2}{3}k_{c}(\unicode[STIX]{x1D6FC})\sim \frac{2}{3}\frac{1}{\cos \unicode[STIX]{x1D6FC}}\frac{{\mathcal{A}}_{0}}{{\mathcal{B}}_{0}}=\frac{1}{\cos \unicode[STIX]{x1D6FC}}k_{0}.\end{eqnarray}$$

Observing that the most unstable mode maximises the growth rate associated with each of the two winds, we obtain

(3.8) $$\begin{eqnarray}k_{max}=k_{\unicode[STIX]{x1D6FC}_{1}(\unicode[STIX]{x1D6FC}_{max})}=k_{\unicode[STIX]{x1D6FC}_{2}(\unicode[STIX]{x1D6FC}_{max})}\end{eqnarray}$$

implying that $\cos \unicode[STIX]{x1D6FC}_{1}(\unicode[STIX]{x1D6FC}_{max})=\cos \unicode[STIX]{x1D6FC}_{2}(\unicode[STIX]{x1D6FC}_{max})$ . The corresponding wavelength is

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{max}\sim \cos \unicode[STIX]{x1D6FC}_{1,2}(\unicode[STIX]{x1D6FC}_{max})\unicode[STIX]{x1D706}_{0}.\end{eqnarray}$$

In particular, for a symmetrical wind regime ( $N=1$ ) and using (3.4), we have

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \cos \unicode[STIX]{x1D6FC}_{1,2}(\unicode[STIX]{x1D6FC}_{max})=\cos \frac{\unicode[STIX]{x1D703}}{2}\quad \text{for }\unicode[STIX]{x1D703}\leqslant 90^{\circ }, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \cos \unicode[STIX]{x1D6FC}_{1,2}(\unicode[STIX]{x1D6FC}_{max})=\sin \frac{\unicode[STIX]{x1D703}}{2}\quad \text{for }\unicode[STIX]{x1D703}\geqslant 90^{\circ }, & \displaystyle\end{eqnarray}$$

which gives

(3.12a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{max}=\unicode[STIX]{x1D706}_{0}\cos \frac{\unicode[STIX]{x1D703}}{2}\quad \text{or}\quad k_{max}=k_{0}/\cos \frac{\unicode[STIX]{x1D703}}{2}\quad \text{for }\unicode[STIX]{x1D703}\leqslant 90^{\circ },\end{eqnarray}$$

(3.13a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{max}=\unicode[STIX]{x1D706}_{0}\sin \frac{\unicode[STIX]{x1D703}}{2}\quad \text{or}\quad k_{max}=k_{0}/\sin \frac{\unicode[STIX]{x1D703}}{2}\quad \text{for }\unicode[STIX]{x1D703}\geqslant 90^{\circ }.\end{eqnarray}$$

For $N>1$ , we do not have an analytic expression for $\cos \unicode[STIX]{x1D6FC}_{1,2}(\unicode[STIX]{x1D6FC}_{max})$ as a function of $\unicode[STIX]{x1D703}$ , but this can be done numerically. These asymptotic predictions of the wavenumber are shown for different values of the transport ratio $N$ in figure 7(a) in dotted lines. We recover with this geometric argument the trends predicted by the linear stability analysis at finite $u_{\ast }/u_{th}$ , and in particular the flattening of the curve for larger values of  $N$ . We notice, however, a slight overestimation of the peak at $\unicode[STIX]{x1D703}=90^{\circ }$ , because the slope effects are still present.

The role of the velocity ratio $u_{\ast }/u_{th}$ on the selected wavenumber is shown in figure 7(b) for $N=1$ . When the wind strength decreases, the threshold and direction of transport are changed as the slope effects are increasingly important (equation (2.9)). First, because the amplitude of the hydrodynamic shift is reduced (the ‘ ${\mathcal{B}}$ effect’), the selected wavenumber is smaller, as in the unidirectional case (2.23). Moreover, for a bidirectional flow regime, $k_{max}$ is even smaller near the transition angle $\unicode[STIX]{x1D703}=90^{\circ }$ . This is a consequence of the stabilising action of the transverse component of sediment transport, as discussed in § 2.2. This effect is the strongest at the transition angle, as it depends on the component of the slope perpendicular to the flow direction. As a result, the curves close to and far from the transport threshold show opposite behaviours. For intermediate wind strengths, the curve $k_{max}(\unicode[STIX]{x1D703})$ is almost flat, as the geometric and slope mechanisms compensate each other.

Figure 8.

Propagation velocity of the most unstable mode $\boldsymbol{c}_{\unicode[STIX]{x1D6F4}}(\unicode[STIX]{x1D6FC}_{max},k_{max})$ and resultant flux $\boldsymbol{Q}_{\unicode[STIX]{x1D6F4}}$ as functions of the divergence angle $\unicode[STIX]{x1D703}$ for $u_{\ast }/u_{th}=1.6$ and various transport ratios $N$ . From dark to light greys: increasing values of $N\in \{1,2,3,+\infty \}$ . (a) Solid lines: propagation velocity $c_{max}$ of the most unstable mode. Dashed lines: resultant flux $Q_{\unicode[STIX]{x1D6F4}}$ . Both are rescaled by their value at $\unicode[STIX]{x1D703}=0^{\circ }$ . (b) Solid lines: direction of the resultant flux $\unicode[STIX]{x1D6FD}_{max}$ measured with respect to the $x$ -axis. For comparison, the pattern orientation $\unicode[STIX]{x1D6FC}_{max}$ (see figure 6 a) is also plotted (dashed lines).

3.5 Propagation velocity of the most unstable mode

Similarly to the growth rate (3.3), the propagation velocity of the emerging pattern can be expressed as the sum of the contribution of each wind:

(3.14) $$\begin{eqnarray}\boldsymbol{c}_{\unicode[STIX]{x1D6F4}}=\frac{1}{N+1}[Nc(Q_{1},\unicode[STIX]{x1D6FC}_{1},k)+\text{sgn}(\unicode[STIX]{x1D703}-90)c(Q_{2},\unicode[STIX]{x1D6FC}_{2},k)]\frac{\boldsymbol{k}}{k}.\end{eqnarray}$$

When the divergence angle is larger than $90^{\circ }$ , the displacements of the pattern induced by each of the winds are opposite, resulting in the term $\text{sgn}(\unicode[STIX]{x1D703}-90)$ . From the propagation velocity can be computed that of the most unstable mode $c_{max}=\Vert \boldsymbol{c}_{\unicode[STIX]{x1D6F4}}(\unicode[STIX]{x1D6FC}_{max},k_{max})\Vert$ . This quantity, as well as the resultant flux $Q_{\unicode[STIX]{x1D6F4}}$ on a flat bed (3.1) and its direction quantified by $\unicode[STIX]{x1D6FD}_{max}=\angle (\boldsymbol{Q}_{\unicode[STIX]{x1D6F4}},\boldsymbol{e}_{x})$ , are shown in figure 8 as functions of the divergence angle $\unicode[STIX]{x1D703}$ for a velocity ratio $u_{\ast }/u_{th}=1.6$ and different values of the transport ratio $N$ .

The propagation velocity is controlled by the component of $Q_{\unicode[STIX]{x1D6F4}}$ perpendicular to the crest and by the pattern wavelength. For $\unicode[STIX]{x1D703}<70^{\circ }$ , the overall behaviour is simple: $c_{max}$ and $Q_{\unicode[STIX]{x1D6F4}}$ vary in the same way. Their directions are almost the same, as the flux direction is not significantly modified by the pattern. As shown in figure 8(a) by the dashed lines, the resultant flux and thus the propagation velocity decrease for wider divergence angles, because the winds increasingly compensate each other. This compensation is weaker when one wind is stronger than the other ( $N>1$ ). For $\unicode[STIX]{x1D703}>110^{\circ }$ , the resultant sand flux is mostly aligned with the dune crest, leading to a slower propagation. The velocity even vanishes for $N=1$ . Finally, around the transition angle $\unicode[STIX]{x1D703}=90^{\circ }$ , the resultant flux direction is strongly modified by the pattern orientation. In addition, the wavelength is also affected by the slope effects. Both effects lead to a complex behaviour of the propagation velocity, which depends on the velocity ratio $u_{\ast }/u_{th}$ .

4.1 Experimental set-up and calibration

The set-up is the one used by Reffet et al. (Reference Reffet, Courrech du Pont, Hersen and Douady2010) and Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014), which we briefly recall here. The granular material (ceramic beads of diameter $d=90~\unicode[STIX]{x03BC}\text{m}$ and density $\unicode[STIX]{x1D70C}=3800~\text{kg}~\text{m}^{-3}$ ) lies on a $80~\text{cm}$ wide square plate, which can move from one side to the other of a $2~\text{m}$ long and $1~\text{m}$ wide water tank. The tank is filled to a water height of $30~\text{cm}$ and the plate stands at mid-height. The plate first quickly moves in one direction at a velocity that generates sediment transport, and slowly translates in the opposite way at $2~\text{cm}~\text{s}^{-1}$ to prevent grain motion. This sequence, which we call a blow, simulates a unidirectional flow over the bed. At the centre of the plate, a $70~\text{cm}$ diameter disc can rotate to change the relative direction of the flow with respect to the bed. A camera set above the plate allows us to follow the evolution of the bed after each blow.

The quick motion of the plate is decomposed into two steps. The plate accelerates from rest to the final velocity in $250~\text{ms}$ and then gently decelerates in approximately $4~\text{s}$ to prevent a backflow. For simplicity, we choose the mean plate speed $V$ as the parameter representing the plate movement. Because the translation of the plate puts the fluid into motion inside a boundary layer that grows over time, the shear stress on the bed is maximum during the acceleration stage and most of the transport takes place during this first phase. The time needed for the boundary layer to reach the water surface is much larger than the translation time, so that the flow can be considered as unbounded.

The transport threshold was not quantitatively determined for this experimental set-up in previous studies. Because our predictions depend on the velocity ratio $u_{\ast }/u_{th}$ , it is important to know how far from the threshold experiments are run. We can provide a first estimation with the observation that no pattern is visible after $2020~\text{blows}$ for a mean plate velocity $V=13.7~\text{cm}~\text{s}^{-1}$ . However, since growth rates are much smaller close to the transport threshold, this estimation is made more precise by using the migration velocity $C$ of barchan dunes of different sizes submitted to a unidirectional flow regime. We vary the mean plate translation velocity $V$ while keeping the plate acceleration time constant. Figure 9(a) shows the velocity of barchans as a function of their width for a fixed translation motion ( $V=20.8~\text{cm}~\text{s}^{-1}$ ). The barchan velocity is smaller for larger barchans, and is well fitted by the Bagnold-like law (Bagnold Reference Bagnold1941)

(4.1) $$\begin{eqnarray}C=\frac{Q_{exp}}{W+W_{0}}.\end{eqnarray}$$

This propagation law, usually expressed in terms of the dune height $H$ , derives from the conservation of mass assuming a constant aspect ratio for the dunes. Under this assumption, the height, the spanwise width $W$ and the streamwise length of barchans are proportional to each other, and $W$ is a more convenient observable on the photos. In (4.1) $W_{0}$ accounts for a correction of Bagnold’s law at small sizes (Elbelrhiti et al. Reference Elbelrhiti, Andreotti and Claudin2008). As illustrated in figure 9(a), this law agrees well with the velocity measurements, and its fit allows us to extract an effective flux $Q_{exp}$ from the different datasets corresponding to various values of $V$ . In this fitting procedure, the adjustable parameter $W_{0}$ has been constrained to be the same for all the data.

Plotting the effective flux $Q_{exp}$ as a function of the plate velocity (figure 9 b), we see as expected an increasing behaviour with a threshold below which the barchans do not move. Although the sediment transport is intermittent in this experiment, we fit the data with the usual transport law $Q_{exp}\propto (V^{2}-V_{th}^{2})$ . This gives a threshold plate velocity $V_{th}=14.2\pm 0.1~\text{cm}~\text{s}^{-1}$ . Note that this value does not depend much on the precise choice of the fitting function.

Figure 10.

Pattern formation from a flat bed in a bidirectional flow regime for $\unicode[STIX]{x1D703}=30^{\circ }$ (left) and $\unicode[STIX]{x1D703}=150^{\circ }$ (right), for $V/V_{th}=1.4$ . (a,b) Pictures after 17 and 13 blows respectively, starting from a flat bed. (c,d) Autocorrelation of the above pictures. Insets: central zone of the autocorrelation matrix encoded in grey scales, which gives the pattern orientation and wavelength. The black curves are the autocorrelation signal along the white dashed lines, i.e. perpendicular to the dune orientation. The measured wavelength $\unicode[STIX]{x1D706}$ corresponds to the position of the first maximum of the autocorrelation curve. (e,f) Wavelength as a function of time. Dotted lines: initial wavelength and corresponding time at which the spatial coherence of the autocorrelation signal extends over one wavelength. Error bars encode the quality of this signal.

4.2 Experimental results

We present here results in the case of symmetric winds ( $N=1$ ): the central disc of the plate is rotated by $\unicode[STIX]{x1D703}$ and $-\unicode[STIX]{x1D703}$ alternately after every blow. Starting from an initial flat bed, dunes rapidly form after a few blows with a characteristic wavelength and orientation, which are measured by picture autocorrelation (see figure 10). We mainly focus below on the dependence of the wavelength with the divergence angle $\unicode[STIX]{x1D703}$ , as the behaviour of the pattern orientation in this experimental set-up has already been discussed in Courrech du Pont et al. (Reference Courrech du Pont, Narteau and Gao2014), showing agreement with their approach (see § 3.3), and thus the present analysis (figure 6).

Figure 10 displays results of two well-contrasted cases with a small and a large divergence angle ( $\unicode[STIX]{x1D703}=30^{\circ }$ and $\unicode[STIX]{x1D703}=150^{\circ }$ , respectively), keeping the same plate velocity ratio $V/V_{th}=1.4$ . When $\unicode[STIX]{x1D703}=30^{\circ }$ , transverse dunes form perpendicular to the resultant transport direction with an initial wavelength of $11.8\pm 0.7~\text{mm}$ . When $\unicode[STIX]{x1D703}=150^{\circ }$ , longitudinal dunes form parallel to the resultant transport direction with an initial wavelength of $8.8\pm 0.5~\text{mm}$ . After the emergence of the pattern, some coarsening is observed, with an increase of the wavelength with time.

Figure 11 shows the behaviour of the wavelength when the divergence angle $\unicode[STIX]{x1D703}$ is varied over the whole range. When $\unicode[STIX]{x1D703}<90^{\circ }$ , transverse dunes form. Their wavelength is maximum for $\unicode[STIX]{x1D703}=0^{\circ }$ and decreases when $\unicode[STIX]{x1D703}$ increases. When $\unicode[STIX]{x1D703}>90^{\circ }$ , longitudinal dunes form and exhibit opposite behaviour: their wavelength increases with $\unicode[STIX]{x1D703}$ and is maximum for $\unicode[STIX]{x1D703}=180^{\circ }$ . Surprisingly, this variation of the dune wavelength with the divergence angle is not symmetric with respect to $\unicode[STIX]{x1D703}=90^{\circ }$ , and longitudinal dunes systematically form at smaller wavelengths than transverse ones. This is true for all values of the plate velocity ratio $V/V_{th}$ that we have tested (see figure 11 a). Furthermore, the wavelength under a unidirectional wind $\unicode[STIX]{x1D706}_{0}$ is found to decrease with $V/V_{th}$ . Wavelengths are significantly larger close to the transport threshold in comparison to the asymptotic value of approximately $12~\text{mm}$ corresponding to large $V/V_{th}$ . By contrast, the wavelength $\unicode[STIX]{x1D706}_{180}$ under two opposed flows ( $\unicode[STIX]{x1D703}=180^{\circ }$ ) is approximately $9~\text{mm}$ and does not depend on the plate velocity.

We can also estimate the characteristic time needed to observe incipient dunes in the experiment. This time is a proxy for the inverse of the instability growth rate. It is displayed in figure 11(b) as a function of the divergence angle $\unicode[STIX]{x1D703}$ . Both transverse and longitudinal dunes emerge typically after a dozen blows, but systematically more rapidly in the case of longitudinal ones. The pattern coarsening is also quicker for longitudinal dunes than for transverse dunes (figure 10(e,f)), an observation already reported by Reffet et al. (Reference Reffet, Courrech du Pont, Hersen and Douady2010). At the transition between the two dune orientations, close to $\unicode[STIX]{x1D703}=90^{\circ }$ , the two modes are observed simultaneously in the building dune pattern, as shown in figure 12. The longitudinal orientation seems to overcome the transverse one at the transition. This asymmetry appears to be specific to incipient dunes, as the transition is sharper once dunes are more mature (Reffet et al. Reference Reffet, Courrech du Pont, Hersen and Douady2010).

Figure 11.

Initial wavelength as a function of the divergence angle $\unicode[STIX]{x1D703}$ (symmetric winds $N=1$ ). Symbols: experimental data for $V/V_{th}=1.4$ . Solid line: model prediction for $u_{\ast }/u_{th}=2.5$ and $L_{sat}=1.1~\text{mm}$ . Dashed line: solid line shifted below and adjusted for $\unicode[STIX]{x1D703}=180^{\circ }$ (we subtracted $(\unicode[STIX]{x1D706}_{0}-\unicode[STIX]{x1D706}_{180})$ ). Inset (a): initial wavelength as a function of the rescaled plate velocity for $\unicode[STIX]{x1D703}=0^{\circ }$ ( $\unicode[STIX]{x1D706}_{0}$ , bullets) and $\unicode[STIX]{x1D703}=180^{\circ }$ ( $\unicode[STIX]{x1D706}_{180}$ , squares). Solid line: fit of the linear stability analysis under a unidirectional flow. It gives $L_{sat}=1.1\pm 0.1$  mm and $C_{v}=4\pm 2$ , see (4.2). Inset (b): characteristic time (in blows) of pattern emergence as a function of $\unicode[STIX]{x1D703}$ . Symbols: experimental data for $V/V_{th}=1.4$ . Solid line: model prediction for $u_{\ast }/u_{th}=2.5$ , adjusted for $\unicode[STIX]{x1D703}=0^{\circ }$ . For all model predictions, ${\mathcal{A}}_{0}=4$ and ${\mathcal{B}}_{0}=4$ (see text).

4.3 Comparison with the theoretical predictions

Figure 12.

Photographs of the dune pattern in the experiment across the transition from transverse to longitudinal orientations, starting from a flat bed, for $V/V_{th}=1.4$ . Flow directions are shown with arrows. From (a) to (c), the pictures have been taken after 30, 32 and 30 blows, respectively.

In order to make comparisons between these experiments and the linear stability analysis developed for aeolian dunes, the model should be adapted to the subaqueous case. The transport law for bedload is different from that of saltation, with a power $3/2$ rather than the linear relationship (2.8) between the saturated sediment flux and the basal shear stress. However, the details of this transport law affect the linear stability analysis by modifying factors in (2.14)–(2.17) only, and we have checked that all the discussed results (e.g. increase of wavelength close to transport threshold, decrease of wavelength close to $\unicode[STIX]{x1D703}=90^{\circ }$ ) are qualitatively robust with respect to any choice of transport law. Also, the transient nature of the flow and grain transport in the experiment prevents us from considering this configuration as a standard bedload situation, where we could quantitatively rely on well-calibrated transport laws. For these reasons, we opted for simplicity and chose to keep the linear transport law (2.8).

Even if we do not necessarily aim at a quantitative matching of the model with the experimental data, the values for the basal shear stress coefficients can nevertheless be more consistent with the experimental case. In fact, the coefficients ${\mathcal{A}}_{0}=3.5$ and ${\mathcal{B}}_{0}=2$ are typical values measured in the field for aeolian dunes of 10 m wavelength (Claudin et al. Reference Claudin, Wiggs and Andreotti2013). These coefficients are weak (logarithmic) functions of the pattern wavenumber at small $kz_{0}$ , but show stronger variations in the range $10^{-4}<kz_{0}<10^{-2}$ due to an anomalous hydrodynamic response of the flow to the bed perturbation, as first experimentally evidenced by Zilker, Cook & Hanratty (Reference Zilker, Cook and Hanratty1977) and Frederick & Hanratty (Reference Frederick and Hanratty1988). In the experiment, the pattern wavelength is of the order of 10 mm, i.e. wavenumber $k=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}\simeq 5.7\times 10^{2}~\text{m}^{-1}$ . The hydrodynamic roughness $z_{0}$ is more difficult to estimate under these transient water flow and sediment transport conditions. Taking it as a small fraction of the grain diameter, we typically have $kz_{0}\simeq 10^{-3}$ , i.e. two orders of magnitude larger than in the aeolian case. For such a value, the coefficients become ${\mathcal{A}}_{0}=4$ and ${\mathcal{B}}_{0}=4$ (Charru et al. Reference Charru, Andreotti and Claudin2013), and we take these values to be fixed for all theoretical computations in this section. Note that increasing these coefficients corresponds to a larger relative weight of hydrodynamic effects compared to slope effects. Interestingly, with these slightly larger values of ${\mathcal{A}}_{0}$ and ${\mathcal{B}}_{0}$ compared to the aeolian case, the slope effects are no longer strong enough to induce an increase of the wavelength at $\unicode[STIX]{x1D703}=90^{\circ }$ close to the transport threshold (figure 7 b). Only the hydrodynamic effect leading to the decrease of the wavelength at the transition angle remains. This is consistent with the experimental observations.

The most unstable wavelength $\unicode[STIX]{x1D706}_{max}$ selected by the linear stability analysis depends on the velocity ratio $u_{\ast }/u_{th}$ . This ratio is not directly equal to the plate velocity ratio $V/V_{th}$ in the experiment because of the transient nature of the flow. Here, we assume that they are linearly related to each other:

(4.2) $$\begin{eqnarray}\left(\frac{V}{V_{th}}-1\right)=C_{v}\left(\frac{u_{\ast }}{u_{th}}-1\right),\end{eqnarray}$$

where $C_{v}$ is an adjustable constant. In the fitting procedure, the most unstable wavelength is thus a function of $V/V_{th}$ with parameters $C_{v}$ and $L_{sat}$ . For a unidirectional flow ( $\unicode[STIX]{x1D703}=0^{\circ }$ ), the comparison between data and theory is shown in figure 11(a). The observed increase of the wavelength close to the transport threshold is quantitatively reproduced by the model. This effect is due to the sensitivity of the threshold to the bed slope. The fit gives $L_{sat}=1.1\pm 0.1$  mm (i.e. $\simeq 12$ grain diameters) and $C_{v}=4\pm 2$ . These values are now fixed and used for the subsequent computation of theoretical curves at all divergence angles, without any further adjustment. In particular, we run all experiments at $V/V_{th}=1.4$ , which corresponds to a velocity ratio $u_{\ast }/u_{th}=2.5$ , from (4.2).

Figure 11 shows that the linear stability analysis captures well the variation of the wavelength for divergence angles smaller than $90^{\circ }$ (solid line). For divergence angles larger than $90^{\circ }$ , the experimental data are significantly smaller than predicted whereas the model predicts a symmetric variation of the wavelength with respect to the transition angle $\unicode[STIX]{x1D703}=90^{\circ }$ . Nevertheless, the decrease of the wavelength close to $\unicode[STIX]{x1D703}=90^{\circ }$ is recovered relative to the wavelength at $\unicode[STIX]{x1D703}=180^{\circ }$ (dashed line). For this configuration of opposite winds, the wavelength is found to be remarkably independent of $V/V_{th}$ (figure 11 a), as opposed to the increase close to the transport threshold observed for $\unicode[STIX]{x1D703}=0^{\circ }$ . Finally, the variation of the characteristic time for pattern formation close to the transition angle is also recovered by the predicted variation of $1/\unicode[STIX]{x1D70E}_{max}$ for divergence angles smaller than $90^{\circ }$ (figure 11 b). However, this is not the case for $\unicode[STIX]{x1D703}>90^{\circ }$ .

These results question the ability to observe the linear regime in the experiment, especially when $\unicode[STIX]{x1D703}>90^{\circ }$ . Such a regime is theoretically associated with bedforms whose wavelength remains constant while their amplitude grows at a rate $\unicode[STIX]{x1D70E}$ . Clearly, as soon as a pattern wavelength is detected by autocorrelation of the experimental pictures, it keeps increasing with time, especially for longitudinal dunes (figure 10 e,f). Enhanced nonlinear effects when $\unicode[STIX]{x1D703}>90^{\circ }$ can be understood as follows. In a bidirectional flow regime, transverse dunes are blown from the same side by the two winds. On the contrary, each wind blows from a different side of the longitudinal dunes. As a consequence, the sand flux is all the more affected if the incipient dunes develop an avalanche face. A larger dune aspect ratio is a stronger perturbation for the flow, leading to larger sediment flux variations, and thus an enhanced dune growth rate. These nonlinearities may explain the shortcomings of the linear analysis in reproducing the experimental wavelengths for $\unicode[STIX]{x1D703}>90^{\circ }$ . These effects also support the shorter characteristic time for dune formation in the longitudinal mode and the prevalence of longitudinal dunes over transverse ones at the transition, when the divergence angle is close to $90^{\circ }$ .

Finally, we can provide estimations of time scales in the experiment, taking as a typical example the case $V/V_{th}=1.4$ . Each blow takes $\simeq 4$  s. Because most of the grain motion occurs during the acceleration phase of the plate, the time during which the dune pattern is active is a fraction of the blow duration. The effective wind period $t_{1}+t_{2}$ is thus of the order of a few seconds. For $V/V_{th}=1.4$ , we have $Q_{exp}=2~\text{mm}^{2}~\text{s}^{-1}$ (figure 9 b). Importantly, this is the flux at the crest of barchan dunes. Taking $\unicode[STIX]{x1D6FE}=1.7$ (see § 3.3), we then obtain a reference flux $Q\simeq 0.74~\text{mm}^{2}~\text{s}^{-1}$ . With a saturation length $L_{sat}=1.1~\text{mm}$ , we have a characteristic pattern growth time $(1/\unicode[STIX]{x1D70E}_{max})L_{sat}^{2}/Q\simeq 5$  s, where $\unicode[STIX]{x1D70E}_{max}\simeq 0.35$ is the dimensionless growth rate given by the linear stability analysis corresponding to the experimental conditions ( ${\mathcal{A}}_{0}=4$ , ${\mathcal{B}}_{0}=4$ and $u_{\ast }/u_{th}=2.5$ ). Despite both time scales being comparable in the experiment, the emerging dunes still develop by integrating the two flow directions. While the experimental conditions are at the limit of validity of the linear stability analysis, some of the characteristics of the linear regime remain, such as the variations of the dune orientation or wavelength with respect to the divergence angle.