Numerical modeling of the wind flow over a transverse dune

Transverse dunes, which form under unidirectional winds and have fixed profile in the direction perpendicular to the wind, occur on all celestial objects of our solar system where dunes have been detected. Here we perform a numerical study of the average turbulent wind flow over a transverse dune by means of computational fluid dynamics simulations. We find that the length of the zone of recirculating flow at the dune lee — the separation bubble — displays a surprisingly strong dependence on the wind shear velocity, u*: it is nearly independent of u* for shear velocities within the range between 0.2 m/s and 0.8 m/s but increases linearly with u* for larger shear velocities. Our calculations show that transport in the direction opposite to dune migration within the separation bubble can be sustained if u* is larger than approximately 0.39 m/s, whereas a larger value of u* (about 0.49 m/s) is required to initiate this reverse transport.

different for isolated or closely spaced dunes 15,16 . However, the effect of the wind speed on the separation streamline is still uncertain. On the one hand, information about the wind shear velocity, u * , is encoded in the shape of the transverse dune. Modeling has shown that the lower u * , the rounder the dune crest 17 , whereas it is known that the zone of recirculating flow becomes smaller the rounder the dune crest 18,19 . On the other hand, the fit parameters for determining s(x) for different dune profiles are typically obtained from simulations performed with one single value of wind shear velocity regardless of the dune shape.
In the present work, we use Computational Fluid Dynamics simulations in order to calculate the average turbulent wind field over a transverse dune of fixed profile under different flow speeds. Specifically, we consider an isolated transverse dune of fixed profile and investigate how the shape and the size of the separation bubble, as well as the characteristics of the flow inside it, vary with the average wind shear velocity, u * , to which the dune is subjected. Figure 2 shows the schematic representation of the setup employed in our calculations. The dune, placed on top of the bottom wall of a two-dimensional channel of height D z 5 20 m and width D x 5 100 m, has height H 5 3.26 m at the brink and width L 5 25 m in the direction of the wind (x). The dune is placed in the center of the box, and the dimensions of the box are chosen in such a way that the dune is far enough from the channel's walls, thus avoiding that the results are affected by border effects. That is, increasing the size of the box does not change the results of our calculations. The wind velocity u(z) increases logarithmically with the height z above the bed level (h) 1,4 . That is, where k 5 0.4 is the von Kármán constant, z 0 is the surface roughness and u * is the upwind shear velocity of the wind, which is used to define the (upwind) shear stress, with r f 5 1.225 kg/m 3 standing for air density. In the present study, we use a value of z 0 5 100 mm. This value has been obtained in a previous work 20 by fitting Eq. (1) to the steady-state wind profile within the simulation box (Fig. 2), which has been generated upon imposing a pressure difference between the extremities of the ''wind tunnel'' using different values of induced flow speed 20 . We find that using other values of z 0 within the range between 10 mm and 1.0 mm 4,21 has a negligible effect on the values of shear stress obtained in our calculations. While previous works investigated the flow over triangular obstacles 11 or artificial dune profiles with a windward side built using halfcircles 16 , the dune profile used in our calculations is generated using a morphodynamic model for dune formation 5,6,14,22,23 . Using this model, we obtain the longitudinal profile of a transverse dune by starting with a Gaussian hill subjected to a wind shear velocity u * 5 0.4 m/s, which is a typical value for the average shear velocity of top. In the bottom-left-hand-corner we see the sand rose for Doha AP in Qatar, located about 170 km SE of the transverse dune field. The vector length of each direction (i) of the sand rose gives the potential rate of sand transport due to sand-moving winds blowing from that direction, i.e. the drift potential, 24,35 . In this equation, U is the wind velocity, U t is the threshold wind speed for transport and f U is the fraction of time the wind speed exceeds U t . The arrow indicates the resultant drift potential, RDP 5 SDP i , where the summation is over all directions of the sand rose.   sand-moving winds on Earth's dune fields 24 . The final dune shape (shown in Fig. 2) has a sharp brink which separates the windward side, on which wind-driven particle transport -that is, saltation 1,25 -occurs, from the slip-face, which has the inclination of the angle of repose 34u and is the side of the dune where avalanches occur (see Fig. 1e).
Motivated by the observation that some dunes on Mars might evolve from sand deposition in the wake of an indurated dune 19 , here we aim to understand how the size of the separation bubble of a transverse dune of fixed profile changes when the wind speed is varied. Therefore, in the present study the same dune profile of Fig. 2 is used to perform all simulations, regardless of the value of the wind velocity.
As mentioned previously, one further aim of the present study is to understand how the flow characteristics inside the separation bubble of a transverse dune changes upon a variation in wind speed. We note that the steady-state profile of a mobile transverse dune is adapting to the flow conditions within the time-scale of dune motion. A systematic investigation of the role of the average wind speed for the longterm evolution of the dune shape would thus require computational fluid dynamic calculations coupled with a morphodynamic model which accounts for the changes in the dune profile due to the winddriven saltation flux. On the other hand, the time-scale associated with the reconstitution of the dune shape after a change in flow conditions -that is, the dune turnover time, T m -is of the order of months or years 4,21,26 , and is thus much larger than the time-scale of the turbulent fluctuations of the wind flow over the dune profile (of the order of seconds 27 ). Thus, keeping the dune profile fixed in the simulations provides an approximation for investigating how changes in the upwind shear velocity occurring at time-scales much smaller than the dune turnover time affect the nature of the flow inside the dune separation bubble 12,28 .

Results
Length of the separating streamline. We calculate the size of the flow separation zone for different values of u * within the range 0.01 m/s , u * , 1.6 m/s. It is useful to express this range of shear velocities in terms of the Reynolds number (Re). In order to define Re, we take the dune height as the characteristic length-scale associated with the flow over the transverse dune 29 . Therefore, we write the Reynolds number as Re 5 U 1 Hr f /g, where H < 3.26 m for the transverse dune investigated here, while g < 1.78 3 10 25 kg m 21 s 21 is the dynamic viscosity of the air, and U 1 is a characteristic value of the wind velocity -which is typically taken at a height of 1 m above the bed 24 . Using Eq. (1) to calculate U 1 5 u(z 5 1 m), we obtain that the above range of u * corresponds to Reynolds numbers in the interval 5.2 3 10 4 , Re , 8.3 3 10 6 . Therefore, all values of u * used in our simulations are associated with fully developed turbulent wind flows, since Reynolds numbers larger than about 6000 lead to turbulent flows in the atmospheric boundary layer 30 . Figure 3 shows the contour plot of the mass flux in the wake of the dune for different values of u * . For all values of u * within the aforementioned range, a well-developed zone of recirculating flow is observed downwind of the dune. This zone, also called ''separation bubble'', separates recirculation from the upper region of laminar flow. The streamline which separates both zones, i.e. the separating streamline, connects the brink to the ground at the reattachment point. The reattachment point is located at position x r , which has a horizontal distance l from the brink position (x b ). We investigate how the size of the separation zone and the shape of separating streamline change with u * .
The position x r of the reattachment point is defined as the downwind position at which the horizontal component of the wind velocity at the ground changes sign from upwind (u , 0) to downwind (u . 0). In other words, the separation streamline is the streamline associated with the smallest mass flux (where the mass flux is in units of kg/s and increases from blue to red in Fig. 3) which does not form a loop. In Fig. 3, the streamline associated with each value of u * is represented by the continuous black line connecting the brink to the ground at position x r . We see that the length of flow separation, l 5 x r 2 x b , increases with u * . This is shown more clearly in the main plot of Fig. 4, which displays the separating streamlines corresponding to different values of u * . We also see in the inset of this figure that the streamlines make an angle h r with the horizontal at the reattachment point, which is different from zero. The fact that the separating streamline makes an angle with the bed at the reattachment point was also observed in previous numerical studies 16 , as well as in experiments on subaqueous dunes 15 .
At the brink of the dune, the separating streamline displays a nontrivial behaviour as shown in Fig. 5. As we can see in this figure, the streamline does not separate from the dune profile at the brink, but at a small distance down the dune slip face. We find that the streamline is about two grid spacings apart the brink position, independently of the grid spacing used. The origin of this behavior is that the steep gradients of the fluid flow at the brink make the simulation innacurate at the brink, and since the grid spacing is a variable associated with the simulation rather than the Navier-Stokes equations describing the flow, we conclude that the dip in the separating streamline near the brink shown in Fig. 5 is a numerical artifact. In a previous work 16 , it was shown that, by considering the part of the separating streamline which curves downwards and neglecting the dip near the brink, the separating streamline can be well fitted by ellipses connecting the brink to the reattachment point at the ground.
In the main plot of Fig. 6 we see the dimensionless reattachment length l/H, where H is the dune height at the brink, for different values of u * . We identify three different regimes for the growth of l/H with u * . The first regime (I) corresponds to u * below 0.2 m/s. The second regime (II) is associated with an intermediate range of u * between 0.2 m/s and 0.8 m/s, while values of u * above this range characterize regime III. We see that for values of u * within regime II (intermediate u * ), the increase of l with u * is small. In this regime, l is around six times H, which is roughly the value observed in previous computational fluid dynamics simulations of the wind flow over isolated transverse dunes 11,16 . However, in both regimes I (low u * ) and III (high u * ) the reattachment length increases with u * much faster than it does within regime II.
It is interesting to note that regime II includes the range of u * that is typical for average sand-moving winds on terrestrial dune fields. In order for aeolian sand transport to occur, the average wind speed must be larger than the threshold for sustained sand transport, u t , which for terrestrial dunes is about 0.2 m/s 1 and thus roughly matches the lower limit of the u * interval corresponding to regime II. Moreover, the upper limit of regime II approximately coincides with the threshold for transport of the sand particles through suspended load, which for terrestrial sand transport is about 4u t (see e.g. Ref. 21). Indeed, since average shear velocities associated with saltation transport and dune formation in terrestrial dune fields are typically smaller than < 0.5 m/s 21,22,24 , previous studies of flow separation were performed with u * values not exceeding this ''upper bound'' 13,16 .
The dependence of the streamline's slope h r at the reattachment point on u * is shown in the inset of Fig. 6. We see that h r is a decreasing function of u * which means that the angle at the reattachment point is larger the smaller the length of the separation zone. Our simulations indicate that, for values of u * corresponding to the average shear velocities of Earth's sand-moving winds (between 0.2 m/s and 0.5 m/s), h r is about 34u, which, interestingly, approximately coincides with the angle of repose of sand (the slope of the slip-face).
Flow characteristics within the separation bubble. In addition to characterizing the separation streamline, it is interesting to study the www.nature.com/scientificreports SCIENTIFIC REPORTS | 3 : 2858 | DOI: 10.1038/srep02858   Fig. 6). In the inset we see that the separation streamlines make an angle h r with the horizontal at the reattachment point, which depends on u * as depicted in the inset of Fig. 6. One aspect which is missing in practically all models of dune formation is that, at the lee of the dune, sand can be transported in the direction opposite to the wind due to the recirculating flow in the separation bubble. As a matter of fact, since net transport within the wake of the dune nearly vanishes, the shear stress within the separation bubble in dune models is usually set as zero 5,14,[31][32][33][34] . However, the shape of dunes at the lee side can be affected by sand transport due to the recirculating wind patterns 4,35 . How does the strength of the reverse flow depend on the upwind wind speed? In order to address this question, we calculate the bed shear stress at the ground within the zone of recirculating flow and in particular we study how the magnitude of the shear velocity associated with the reversing flow within the separation bubble depends on the upwind wind shear velocity, u * . The results of our calculations are shown in Fig. 7 and discussed in the paragraphs which follow.
By neglecting the lateral component of the shear stress (that is, the one parallel to the axis of the transverse dune), we can write, whereê x is the unit vector that points in wind direction. In the upper inset of Fig. 7 we show the value of t x within the separation bubble as a function of the horizontal distance x relative to the brink position x b for different upwind shear velocities u * . In the simulations, t x is computed with Eq. (2) by using the local wind shear velocity u * associated with the velocity profile u x (z) (where u x is the component of the wind velocity in the x direction). We see that t x has negative values within the separation bubble. The behavior of t x near the brink is affected by the strong gradients of the flow resulting from the discontinuity of the dune surface at the brink, as discussed previously. Moroever, the behavior of t x inside the separation bubble is sensitive to the strong gradient of h(x) that occurs at the position of the slip face foot, that is, at x 5 x f < 26.3 m. As we can see in the upper inset of Fig. 7, for all values of u * , t x has a dip at that position (that is, at After position x f , t x curves downward until it achieves a minimal value. We define l rev as the distance x 2 x b associated with the smallest value of t x inside the separation bubble, which we call t rev . For instance, in the upper inset of Fig. 7, the point (l rev , t rev ) that corresponds to the upwind shear velocity u * 5 1.0 m/s is denoted by a black dot. Note that while the dip of the separation streamline two grid spacings after the dune brink (cf. Fig. 5) is a numerical artifact which occurs due to the flow separation at the sharp brink but is independent   In the upper inset, we see the value of t x (cf. Eq. (3)) associated with the reversing flow within the separation bubble as a function of the downwind distance x 2 x b , where x b is the brink position. The minimum of t x , which gives the maximal magnitude of the shear stress associated with the reversing flow, is denoted by t rev and occurs at a horizontal distance l rev downwind of the brink. The blue horizontal line indicates the impact threshold shear stress, t t~rf luid u 2 t . The black dot identifies the point {l rev , t rev } associated with u * 5 1.0 m/s. We see in the lower inset that l rev < 55% of the reattachment length l and is nearly independent of u * . Main plot: Magnitude of the maximal shear velocity at the ground, u Ãrev~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi t rev j j=r f luid p , as a function of the upwind shear velocity, u * . The best fit using the equation u *rev 5 ku * (dashed line) gives k < 0.55. The dashed area indicates the range of u * for which no sustained transport occurs (u *rev , u t ), while the regime of u * where u *rev is smaller than the threshold for direct entrainment (u ft ) is denoted by the green filled area.
www.nature.com/scientificreports SCIENTIFIC REPORTS | 3 : 2858 | DOI: 10.1038/srep02858 of the dune shape 16 , the profile of the shear stress is affected by the local slope of the surface. That is, both the dip of the shear stress occurring at the downwind foot of the dune and the saddle of t x at position x b 1 l rev (cf. upper inset of Fig. 7) are no numerical artifacts but result from the local gradients in the slope of the envelope h(x) comprising the dune surface and the ground.
From the magnitude of t rev , we obtain the value of the shear velocity u Ãrev~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi t rev j j=r air p that is associated with the reversing wind flow at the position l rev . The values of u *rev obtained for different values of the upwind shear velocity u * are displayed in the main plot of Fig. 7 (squares). The dashed line denotes the best fit to the simulation data using the equation, which gives k < 0.55. The lower inset of Fig. 7 shows the ratio l rev /l as a function of u * . We see that the maximal value of u *rev associated with the recirculating flow occurs at a downwind distance l rev that is approximately equal to 55% of the reattachment length l, whereas only a weak dependence of l rev /l on u * is observed in our simulations. Based on the results for u *rev found from our simulations, we can estimate for which values of the upwind shear velocity u * surface transport of sand within the separation bubble can occur in the direction opposite to the dune migration trend. In order for sand transport to begin, the local wind shear velocity must exceed the socalled fluid threshold u ft . This threshold shear velocity is given by the equation 1,21 , where g is gravity, d and r p are the mean diameter and density of the sand grains, respectively, while A is the Bagnold-Shields parameter, which depends on the shape and sorting of the grains and on the angle of internal friction 29,36 . Furthermore, the value of A depends on the attributes of sediment and fluid and can be estimated by numerically solving the following empirical equation derived by Iversen for Re *ft $ 10, where Re *ft 5 u ft d/n is the friction Reynolds number associated with the threshold shear velocity u ft , while the constant 6.0 3 10 27 has units of kg m 0.5 s 22 . By considering that the average grain size of the sand that composes terrestrial dunes is d < 250 mm 4 and taking r p < 2650 kg/m 3 (density of quartz), we obtain A < 0.118 and, which gives the minimal value of u *rev required to initiate transport in the direction opposite to dune migration within the separation bubble. Indeed, once sand transport begins, it can be sustained at shear velocities lower than u ft due to the splash mechanism resulting from the grain-bed collisions. The minimal shear velocity required to sustain saltation transport is the so-called impact threshold shear velocity, u t . As found experimentally and in numerical simulations, u t is approximately equal to 80% of u ft 1,21 . Therefore, using the values of d and r p associated with the characteristics of terrestrial dunes' sand specified above, we obtain, This value of shear velocity gives, thus, the minimal value of u *rev required for reverse transport within the separation bubble to be sustained. We find that the values of upwind shear velocity u * associated with the minimal values of u *rev required to initiate and sustain transport in the direction opposite to dune migration within the separation bubble are u *u < 0.49 m/s and u *b < 0.39 m/s, respectively. The dashed area in the main plot of Fig. 7 denotes the range of u * where no transport can occur (u *rev , u t ). Furthermore, the green filled area identifies the regime of u * where transport within the bubble can be sustained, once initiated (u t , u *rev , u ft ), whereas transport initiation requires u * to exceed the upper limit of this intermediate range.

Discussion
This work presents the first systematic calculation of the turbulent wind flow over a slice of a transverse dune for different values of the wind shear velocity, u * . Using computational fluid dynamics, we have shown, for the first time, that the length of the separation streamline at the lee of the dune increases with u * , and that the separation streamline has an angle with the ground at the reattachment point which decreases with u * . These results will have implications for morphodynamic modeling of sand dunes, since the shape of the separating streamline affects the analytical solution for the shear stress over the windward side of the dune 5 .
Moreover, our calculations allowed us to estimate, for the first time, the values of upwind shear velocity required to induce transport of sand in the direction opposite to dune migration within the separation bubble. It is important to remark that our estimation of the minimal upwind shear velocity required to induce sustained saltation within the bubble (that is, our estimation of u *b ) did not account for the transient distance required for the saltation flux to adapt to a change in flow conditions. This transient distance, called saturation length (L sat ), is of the order of 50 cm and depends only marginally on the wind velocity 38 . Therefore, saltation transport in the reverse direction within the separation bubble can be sustained only if the magnitude of t x is above the impact threshold shear stress, t t <r f luid u 2 t , over a distance of about 50 cm. Indeed, we find that, for an upwind shear velocity u * 5 0.4 m/s, the value of t x varies by only 0.8% over a distance of 50 cm upwind from the position l rev . In particular for this value of u * , which is close to our estimate of u *b , jt x j is above t t over an upwind distance of about 3 m from the position l rev , that is almost six times L sat . Therefore, we conclude that including the effect of the saturation length affects only marginally our estimation of the lower bound for the green area in Fig. 7.
By using Eq. (1), we obtain that the wind velocity required for inducing sustained (reverse) transport in the separation bubble is about 5.0 m/s (at a height of 1 m), while an upwind flow velocity of about 6.3 m/s is required for initiating reverse transport in the dune wake. Our results shed light on the occurrence of sand transport gusts in the direction opposite to dune migration trend within the separation bubble, as reported from field observations 12,28 . In fact, the wind velocity can vary considerably over time, as illustrated by the measurements of Ref. 27. These authors reported measurements of the surface wind speed within a barchan dune field over about 2 hours using an anemometer which resolved changes in wind speed within intervals of 1 s. The authors found strong, rapid fluctuations in wind velocity within the broad range between 4 m/s and 11 m/s 27 . Such time fluctuations of the upwind shear velocity combined with the so-called hysteresis effect in sand transport 39  -provide an explanation for the sporadic gusts of sand transport within the separation bubble.
We note that the relevance of wind fluctuations for the occurrence of transport in the separation bubble should depend on the shape of the dune and also on the values of upwind shear velocity, which can further change over different dune fields. Moreover, the geomorphological implications of the reverse transport within the separation bubble are still uncertain 12,28 and should be investigated with modeling. Indeed, in real conditions the wind is not only varying in speed but also in direction 28 . Thus, future work is required for improving our understanding of the secondary flow patterns on the wake of three-dimensional dune shapes, including the effect of wind directions oblique to the symmetry axis of the dune.
It would be interesting to perform the calculations of the present work considering three-dimensional dune shapes and to investigate the role of u * for the patterns of secondary flow in the lee of barchans and transverse dunes. Moreover, the calculations should be repeated for Martian conditions, in particular for investigating the effect of u * on the reversing flow within the separation bubble of Martian dunes. Indeed, the value of u ft /u t on Mars can be much larger than it is on Earth, and can reach values close to ten 39,40 . Future modeling should account for results reported in this paper in order to investigate their implications for the dune morphodynamics.
In particular, understanding some of the exotic dune types occuring in extraterrestrial environments requires an accurate description of the shape of the separation bubble and the wind profile within the zone of recirculating flow 21 . In the north polar region of Mars known as Chasma Boreale, unusual dune forms occur which do not find terrestrial counterparts. The unusual dune morphologies are straight linear dunes and rounded barchans or ''dome-like'' dunes with incipient slip-face, both dune shapes occurring side-by-side (cf. Fig. 8a). From the orientation of the barchan (or dome-like) dunes in Chasma Boreale, it is clear that the straight dunes are longitudinal dunes, that is, they elongate along the wind direction. However, on Earth, longitudinal dunes form in areas of bimodal wind regimes and usually display a characteristic meandering 4,35 . Moreover, it was shown that a longitudinal dune subjected to an unimodal wind is unstable and decays into a longitudinal chain of barchans 26,45 . Thus, if the wind regime at Chasma Boreale is unimodal (which is the wind regime consistent with the occurrence of the barchans), then the linear dunes shouldn't occur.
Recently it was suggested that a different mechanism might be at play in the formation of the dunes at Chasma Boreale, namely sand induration or cementation by salts or moisture in the dune 19 . According to this hypothesis, the straight linear dunes at Chasma Boreale would form from deposition of sand in the wake of the (indurated) barchan dunes. More precisely, according to this picture, sediment transported from the upwind over the windward side of the (non-erodible) barchan is trapped within the separation bubble. As the deposited sediment undergoes cementation, the dune elongates. However, morphodynamic simulations 19 produced indurated domelike dunes which did not elongate further to form straight linear dunes.
While the simulations of Ref. 19 were performed with a constant wind speed, here we tested the hypothesis that an indurated domelike dune could elongate into a straight longitudinal dunes if the wind strength (and thus the separation bubble length) oscillates in time. Specifically, the volume of sand accumulated in the separation bubble during the period of strong wind could serve as source of sand for downwind transport and elongation of the dune in the direction parallel to the wind during the period of low wind strength (and thus short separation bubble length).
However, we could not produce elongated dunes using twodimensional simulations of the dune's longitudinal profile with oscillating wind strength. The initial profile in our simulation is the transverse dune of Fig. 4, which we use to compute the separation streamline for a large value of u * . Next, we take the envelope comprising the dune profile and the separation bubble as a new (indurated) sand surface, for which we then compute the average turbulent wind flow using a much lower value of u * . However, since this new dune profile is smooth (cf. Fig. 4), flow separation does not set in anymore, regardless of further oscillations in the wind speed. The predictions of our simulation is that the final dune shape is the rounded (dome-like) dune (Fig. 8b), rather than the straight dune.
Indeed, it was also conjectured in Ref. 35 that three-dimensional flow patterns -due to deflection of the flow on both flanks of the dune into the separation bubble -at the dune lee contribute an essential factor for driving the elongation of the dune into a straight longitudinal dune. The result of our calculations apparently corroborate the hypothesis of Ref. 35 that the elongation of the Chasma Boreale straight dunes cannot be understood from two-dimensional calculations (in contrast to the formation of transverse dunes). Therefore, future modeling of the Chasma Boreale straight linear dunes should be conducted using calculations of the turbulent wind flow over the three-dimensional dune shape following the present work. we see the transverse dune profile with the height h rescaled by H (the height at the brink), whereas the downwind position x is rescaled by the dune's along-wind width L. In (d) we see the surface which results after the sand incoming from the upwind is deposited within the separation bubble calculated with u * 5 0.9 m/s. By assuming that this surface is indurated, a rounded dune shape is obtained, which does not elongate further to form a straight dune. Methods In our simulations, we consider that the fluid (air) is incompressible and Newtonian, which is a reasonable assumption since the flow velocities involved in our simulations are much smaller than the speed of sound, whereas the average turbulent wind field over the terrain is computed as described in previous works 16,18 . We adopt the FLUENT Inc. commercial package (version 6.1.25) in order to solve the Reynoldsaveraged Navier-Stokes equations with the standard k{E model, which is used to describe turbulence.
Boundary conditions. The logarithmic wind profile given by Eq. (1) is imposed as boundary condition at the inlet of the channel, and the shear velocity u * is the only parameter we change in the calculations. At the outlet of the channel, a constant pressure (P 5 0) is imposed in order to generate a pressure gradient in the flow direction. The no-slip boundary condition is applied to the entire solid-fluid interface comprising the dune surface and the bottom wall, whereas for the top wall we set both components of the shear stress equal to zero, which means that the top wall moves downwind with velocity equal to the flow velocity at the height of the wall (see also Refs. 18,20,40).
Discretisation scheme and specification of the turbulence model. We solve the time-averaged (or Reynolds-averaged) Navier-Stokes equations for the wind flow over the dune in the fully-developed turbulence regime. We use the so-called k{E model with renormalization group (RNG) extensions, which is known to yield most accurate results for flow separation 12,16,41,42 . We choose the default pressure-velocity coupling scheme (''SIMPLE'') of the FLUENT software and use its preselected values of parameters. We also choose the software's default option ''standard wall functions''. These wall functions apply the wall boundary conditions to all solution variables of the k{E model that are consistent with the logarithmic law for the wind velocity (Eq. (1)) along the entire bottom wall of the channel 43 . We employ the second-order upwind scheme for the turbulence kinetic energy, turbulence dissipation rate and velocity, while for the pressure the second-order discretisation scheme is adopted. We use a triangular grid with average spacing 5 cm, which is refined near the dune-fluid interface and in the wake region behind the dune close to the ground.
The initial condition is such that pressure and velocity are set to zero for all values of x and z, except at the left wall (x 5 0), at which the logarithmic profile for the wind velocity, i.e. Eq. (1), is imposed. The transport equations for the RNG k{E model 44 are then numerically solved iteratively until the convergence criteria are fulfilled. In our simulations, the convergence criteria are defined in terms of residuals, which give a measure of the degree to which the conservation equations are satisfied throughout the flow field. We consider that convergence is achieved when the normalized residuals for both velocity components fall below 10 27 and when the normalized residuals for both k and E fall below 10 24 .