Particle dynamics of a cartoon dune

The spatio-temporal evolution of a downsized model for a desert dune is observed experimentally in a narrow water flow channel. A particle tracking method reveals that the migration speed of the model dune is one order of magnitude smaller than that of individual grains. In particular, the erosion rate consists of comparable contributions from creeping (low energy) and saltating (high energy) particles. The saltation flow rate is slightly larger, whereas the number of saltating particles is one order of magnitude lower than that of the creeping ones. The velocity field of the saltating particles is comparable to the velocity field of the driving fluid. It can be observed that the spatial profile of the shear stress reaches its maximum value upstream of the crest, while its minimum lies at the downstream foot of the dune. The particle tracking method reveals that the deposition of entrained particles occurs primarily in the region between these two extrema of the shear stress. Moreover, it is demonstrated that the initial triangular heap evolves to a steady state with constant mass, shape, velocity, and packing fraction after one turnover time has elapsed. Within that time the mean distance between particles initially in contact reaches a value of approximately one quarter of the dune basis length.

" 'We are right in the open desert,' said the doctor. 'Look at that vast reach of sand! What a strange spectacle! What a singular arrangement of nature!' " [1]. In this fictive story of Jules Verne from the 19th century, Dr. Samuel Ferguson and his two companions discover the beauty of sand dunes in the deserts of Africa by traveling five weeks in a balloon across the continent. Today, scientist are still overwhelmed by these self-organized granular structures, but use satellites for observations instead of balloons [2,3].
The more descriptive science of the past century [4,5,6] has changed to a detailed approach in understanding the dune dynamics [7]. The crescent-shaped barchan dune was chosen as a suitable object because of its relatively fast dynamics. So-called 'minimal models' were established to describe the basic dynamics of barchan dunes [8,9,10]. These two-dimensional models deal with dune slices along the direction of the driving wind. Shortly speaking, they combine an analytical description of the turbulent shear flow over low elevations [11,12] with a continuum description, which models the saltation on the surface of the dune [13]. In the next step of complexity, these twodimensional slices are coupled in the cross-wind direction to model three-dimensional barchan dunes [14]. Laboratory experiments [7,15,16,17,18] and field measurements on Earth [19,20,21] or satellite observations of Mars [22,23] reveal the quality of such models. Some aspects of the minimal models have been checked quantitatively in a narrow water flow channel [24,25,26], and the existence of a shape attractor for barchan dunes [27,28] has been demonstrated experimentally [29].
Minimal models are continuum models: They deal with the overall shape of the dunes neglecting the particulate nature of these granular systems. However, the dynamics of dunes is determined by the transport mechanism of the individual grains of sand [4]. The aeolian sand transport consists of two modes of transport: reptation in the low energy regime and saltation in the high energy regime [30]. The saltating grains are carried with the wind and have flight lengths of thousands of grain diameters. The paths of the reptating grains are much shorter. Their motion is initiated by impacts of the saltating grains onto the granular bed [31]. Experiments in wind-tunnels with sand traps [32,33,34] and with particle tracking methods [35,36,37] give insight into the aeolian grain transport on the level of particle dynamics.
Neglecting the possibility of suspension at high shear velocities, the bed-load transport in water can also be separated in two energy regimes: saltation and surface creep [38,39,40,41]. The saltating particles have much shorter flight lengths than in air and, in contrast to the reptating particles, the creeping particles are directly dragged by the fluid and stay always in contact with the bed surface. Particle tracking experiments in flow channels with denser fluids than air, like water or silicon oil, have been performed to investigate the erosion, the transport, and the deposition of grains of sand on a granular bed [42,43,44,45,46].
The particle dynamics at the surface transfers itself into the bed. The packing density of the granular bed decreases from the inner layers towards the outer fluidized layer, which has a thickness of a few grain diameters in a laminar flow [45]. The observation of segregation effects indicate that the sediment in the sand bed is mixed during the process of the sand transport and the associated ripple formation [47]. Besides the global granular transport, the transition between creeping and saltating particles is also a matter of interest [32,42,43,48,41].
We investigate a migrating isolated dune on a plane surface reminiscent of a barchan dune in the desert. Our experiment is supposed to reveal the details of the particle transport above, on, and inside this solitary dune. Together with the characterization of the driving water flow field, the results should encourage future particulate models and simulations [49]. To get access to the dynamics of the individual grains, we create an experimental realization of a minimal dune model -a cartoon dune.
A sketch of our experimental setup is shown in Fig. 1. The main part consists of a closed flow channel machined from perspex. The height of the channel amounts to H = 60 mm, and its width is 50 mm. The length of the straight section is 600 mm, the curves have an outer diameter of 500 mm and an inner diameter of 400 mm. The channel is filled up with distilled water. The flow is generated by a propeller with a diameter of 45 mm, which is installed in the curve following the section of measurements. The propeller is driven by a motor with a shaft which is placed in the middle of the channel profile. The flow direction is counter-clockwise on top view. We narrow the section of measurements to create quasi two-dimensional conditions, comparable to a Hele-Shaw cell [50]. For that purpose we use a black plastic insert, which constricts the channel to a width W = 3 mm. This leads to an aspect ratio of H/W = 20.
For the calculation of the Reynolds number Re in the 3 mm wide gap we measure the flow velocity u with an ultrasonic Doppler velocimeter (UDV) manufactured by Signal Processing SA. This device measures the vertical velocity profiles u(y). From these profiles we extract the mean horizontal velocity u mean by averaging along the central 80 percent of the profiles. The kinematic viscosity of water at experimental temperature (20.5 ± 0.3 • C) is ν ≈ 1 mm 2 /s. The flow velocity is kept at u mean = 0.65 m/s, which corresponds to Re = u mean H/ν = 39000. Instead of natural sand we use spherical white glass beads (SiLibeads), which have a radius of r = 1.00 ± 0.01 mm and a density of ρ = 2.51 g/cm 3 . For the model dune, the overall mass of the beads amounts to m = 6.5 g, which corresponds to 629 particles.
To characterize the particle transport in our experiment, we investigate the threshold of incipient motion for the used glass beads. Therefore, we prepare a flat granular bed and determine the shear velocity near the surface. The corresponding profile u x (y) and the slope within the laminar boundary layer are shown in Fig. 2. We observed that particle motion sets in at the critical shear velocity u * c = 14 mm/s. This number is used for the calculation of the critical particle Reynolds number Re c = u * c 2r ν and the critical Shields parameter with the acceleration of gravity g = 9.81 m/s 2 , the density of water ρ f = 1 g/cm 3 , and the specific density s = ρ/ρ f . The resulting Re c = 28 is about six times larger than the typical value for aeolian dunes Re c, air = 5, whereas Θ c = 0.007 is comparable to its aeolian counterpart Θ c, air = 0.010 [29].
The model dune is observed with a high-speed camera (IDT, MotionScope M3) which is placed in front of the straight part of the channel. The camera is set to a resolution of 1280 × 256 pixels at 2000 frames per second. For a sufficiently bright illumination we use two bars of light-emitting diodes (LED-bars), as shown in Fig. 1. They lighten the section of measurements from below and from above to obtain a homogenous brightness.
The experimental procedure is as follows: After the flow tube is filled with distilled water, a funnel with a 3 mm long slit is used to pour glass beads into the channel. The experiment starts with the triangular heap shown in Fig. 3 at t = 0 s. The five subsequent snapshots show the temporal evolution of the initial triangular heap towards an asymmetric heap moving downstream. After t = 3 s the heap reaches the characteristic steady-state shape known from central slices of barchan dunes along their migration direction [19]. In the images the glass beads appear as bright disks in front of the dark background. For their localization we use the Circle Hough Transformation as described in detail in Refs. [51,52]. With this image-processing technique all bead positions in every recorded frame are found. The high-speed camera allows particle tracking for each individual grain, which makes it possible to study the particle dynamics of dunes and the mixing process within dunes. To visualize the mixing we provide Fig. 4: At the beginning all detected beads of the data set are colored with a continuous color gradient from the left to the right, as the snapshot at t = 0 s shows. As time evolves the traceable beads keep their initial color, the non-traceable ones are painted black. As the flow sets in, the grains at the surface of the upstream-side are carried away and fall onto the downstream-side, where they roll further downwards. In nature this side is called the slipface. From the color code, it can be seen that the inner beads rest in place. The first cycle of the mixing is passed after t = 5 s, when again most of the blue circles are at the windward side.
For a quantitative measure of the mixing process we investigate the temporal evolution of the distances d i,j between those particles i and j which are initially in contact. Six representative histograms are plotted in Fig. 5. The initially sharp distribution broadens with time, and after t = 5 s the distribution spans the complete basis dune length l dune . This basis length is determined by the bulk part of the dune as defined later (see Fig. 12(a)) and turns out to be l dune = 190 mm. The growth of the mean distance d =< d i,j > and its standard deviation is shown in the insert of Fig. 5. It indicates that within the first turnover time t turnover of about 5 s the mean distance between particles initially in contact reaches a value of approximately one quarter of the dune basis length. Moreover, once all particle positions are known, other quantities describing the dynamics of the dune can be extracted. As shown in Fig. 6(a) the number of glass beads N is slightly fluctuating within one percent. This is primarily due to particles flowing in and out of the frame, and partly caused by detecting failures. Although only the particles from the surface are in motion, the center of mass of the dune migrates downstream with a velocity v com . From the positions of all beads in the frame, both the center of mass and its velocity are determined. It can be seen in Fig. 6(b) that v com takes a constant value after t = 2 s. This coincides with the relaxation towards the steady-state shape, which is demonstrated by the constant height of the center of mass h com in Fig. 6(c).
Alternatively, the relaxation process can be characterized by measuring the area packing fraction ϕ a within the core of the dune, which is composed of the inner static particles. For that purpose we calculate the Voronoi tessellation of each frame [53]. The packing fraction is obtained as the ratio between the cross sectional area of a spherical glass bead and the corresponding area of its Voronoi cell. In order to characterize only the core beads, we exclude those beads with ϕ a < 2/3 ϕ max . The maximal packing fraction ϕ max can be larger than one, because an overlap of particles is possible in the two-dimensional projection of a three-dimensional packing. In our narrow channel of thickness 3r the maximal packing fraction ϕ max = π/3 is obtained by assuming a twodimensional square lattice with a lattice constant of 2 √ 3r. The crystal basis consists of two beads in contact forming an angle of 30 • with respect to the plane of this lattice.
In Fig. 7 the resulting Voronoi cells of the six snapshots of Fig. 3 are overlayed to all detected beads drawn as grey disks. The cells are color coded with respect to their area. The loss of the blue cells indicates that the packing fraction decreases during the time of measurements. This is quantitatively demonstrated by Fig. 6(d), where the temporal evolution of the spatial mean packing fraction ϕ mean corresponding to these Voronoi cells is plotted. In addition to the relaxation of v com and h com , ϕ mean also takes a constant value after t = 3 s.
The investigation of the transient demonstrates that the initial triangular heap evolves to a steady state with constant mass, shape, velocity, and packing fraction after about one turnover time t turnover has elapsed. Within that time the mean distance between particles initially in contact reaches a value of approximately one quarter of the dune length l dune .
The presented Voronoi method has the advantage to reveal the temporal evolution of the packing fraction, but has the disadvantage to be operational only within a small fraction near the core of the dune. In order to get meaningful results also in the surrounding area of the dune, we obtain the local area packing fraction from the possibility for each pixel to be shaded by a particle during a certain time. The relation between ϕ a and the volume packing fraction ϕ v is where N is the number of particles in a certain area A. In our geometry r/W = 1/3 and thus ϕ v = 4/9 ϕ a . This yields for the maximum value 4/9 · π/3 ≈ 0.47. The experimentally obtained field ϕ v (x, y) is shown in Fig. 8. Here, the time average is taken over one second, which corresponds to 2000 single frames. For technical reasons the frames are shifted according to the center of mass of the dunes.
The mean velocity field of the particles is calculated similar to ϕ v (x, y). Its absolute value v(x, y), its x-component v x (x, y), and its y-component v y (x, y) are plotted in Fig. 8. The field v(x, y) shows that the velocity increases with increasing height above the dune surface. This indicates that the grain movement changes from creeping to saltationi.e. freely swimming in the water stream. This behavior is also reflected by v x (x, y). Notably, v x (x, y) has almost no negative values, indicating that the recirculation bubble at the downstream side [54] is not strong enough to push the sinking beads backwards. The field v y (x, y) visualizes the erosion and deposition process: at the upstream side the particles are lifted and above the slipface they rapidly fall down.
The motion of the individual particles results from the interplay with the driving water flow. To illuminate this interaction we need the velocity field of the water flow, which can be measured with the UDV. However, the UDV needs a couple of seconds to acquire a single velocity profile. This is too long for recording the velocity field within the time the dune needs to propagate along the section of measurements. We evade this problem by measuring the flow field around a steady dune made of plastic instead of the moving dune. The contour of this dummy is obtained from the snapshot at t = 4 s in Fig. 3 by drawing a smoothed envelope around the densely packed bulk of the dune. The result is shown as a red line in Fig. 9(a). The dummy is decorated with glass beads to simulate the surface roughness, which determines the size of the laminar boundary layer. We define five characteristic areas (A -E) as shown in Fig. 9(a) to compare the vertical velocity profiles of the grains and the water flow in Fig. 9(b). The vertical grain velocity profiles v x (y) are extracted from v x (x, y) by averaging horizontally in the corresponding areas. The vertical water velocity profiles u x (y) are recorded by crossing two ultrasonic beams above the dune surface to get an averaged data set within each area.
The highest measured velocities are found near the crest of the dune. For the water flow this is a manifestation of the incompressibility condition in our closed channel. It is notable that in the areas A to D the particle velocity is finite at such positions y where the water velocity is zero. This can be explained by the fact that the flow velocities are measured with the plastic dummy, which fails to take into account the  . Panel (a) shows the plastic dummy decorated with glass beads and a detail of the snapshot at t = 4 s in Fig. 3. The red line indicates the profile of the plastic dummy used for the determination of the water flow field. The five areas tagged with white lines and capitals correspond to panel (b), which provides the comparison between the vertical velocity profiles v x (y) of grains (red lines) and u x (y) of water flow (blue lines) in the selected areas. Each profile consists of 100 measured data points of which a few representatives are plotted with error bars. motion of the particles on the surface. In area E the water flow velocity changes its sign at lower heights, which is an indication of the recirculation bubble. In contrast, the grain motion does not follow the water flow in the wake of the dune. In all areas the velocity gradients of the water profiles are steeper than the ones of the beads. This might result from the fact, that the particles cannot follow the water flow immediately due to inertia. Notably, the maximal velocities of both, water and beads, is one order of magnitude greater than the migration velocity of the dune measured as the center-of-mass motion. This is explained by the particle size compared to the dune size.
From the velocity gradient of the water flow near the surface of the dune(see Fig. 2), Figure 10. Shear stress of the water flow at the dune surface. A smoothing of the data points yields the red line, which is supposed to guide the eye. The grey shaded area symbolizes the location and the shape of the plastic dummy.
we calculate the shear stress of the water flow τ = η du x dy surface assuming η = 10 −3 Pas for water at 20 • C. Here, the additional turbulent shear stress is neglected, because near fixed walls in the area of the viscous boundary layer the eddy viscosity goes to zero [55]. The resulting values for τ as shown Fig. 10 quantify the slope of the water profiles near the surface of the dune and give a measure for the local drag force. The shear stress increases along the upstream side of the dune and reaches a maximum value afore the crest. Above the slipface τ decreases and even reverses sign in the region of the recirculation bubble and relaxes further downstream. This qualitative behavior corresponds to the analytical description of the wind shear stress used in the theoretical model by Kroy et al [9,10]. The resulting particle fluxes j x (x, y) and j y (x, y) are shown in Fig. 8. They are the product of the particle density ϕ v (x, y)/V bead times v x (x, y) or v y (x, y), respectively, where V bead = 4π/3 mm 3 is the volume of a glass bead. In analogy to Fig. 9(b) we extract the averaged vertical flux profiles of the horizontal flux j x (y) as presented in Fig. 11(a). From these profiles the particle flow rates as a function of the vertical position are extracted by integrating j x (y) with the vertical position y and the width W of the flow channel. The resulting curves are plotted in Fig. 11(b).
To specify the nature of the particle transport, we distinguish between beads which have contacts to other beads forming the bulk shown as gray shaded area in Fig. 12(a), and beads which move freely in the water stream without contact to the bulk. The contact criteria is fulfilled, when the centers of two beads take a distance shorter than  Figure 11. Panel (a) shows the particle flux j x (y) and panel (b) shows the particle flow rate I x (y) for all beads (black lines), beads with contacts (red lines), and freely moving beads (blue lines). The data sets are obtained from the corresponding five areas in Fig. 9(a).
2.1r. The beads within the bulk are composed of resting beads and creeping beads, but only the creeping beads contribute to the flux density. The motion of the free flying particles above the surface is named saltation. It can be seen that at every position along the dune the average vertical positions of creeping particles is near the surface h(x) and the average vertical positions of saltating particles is located two to three particle diameters above. The vertical positions are determined as the first moment of the vertical particle flux distributions. The distributions are calculated with a spatial running average along the dune with a fixed width of 13 mm corresponding to the areas in Fig. 8.
We define the values of the spatially averaged integral particle flow rate I x,creep = I x (y = H) for creeping particles. I x,salt for saltating particles is defined likewise. The total particle flow rate results as I x,total = I x,creep + I x,salt . All three flow rates are shown as a function of the horizontal position in Fig. 12(b). The spatial profiles of the separated flow rates show, that the saltating beads contribute a larger part of the dune motion than the creeping beads. Both flow rates reach their maximum values near the crest of the dune, as it is known from natural dunes in the desert [13,56].
The overall erosion rate q = d 2 N dtdx can be obtained as the spatial derivative of the total flow rate dI x,total /dx. The erosion rate as shown in Fig. 12(c) matches qualitatively with the local slope of the dune surface dh/dx. If q > 0 the beads are entrained into the water stream and will strengthen the total flow rate. If q < 0 the beads are deposited on the surface of the dune and the value of I x,total will decrease. Along the upstream side q is greater than zero and becomes negative behind the crest with a maximal deposition onto the slipface. The plateau near the crest indicates a constant flow rate. These observations correspond to the color code of j y (x, y) in Fig. 8. Notably, the position of the maximal deposition matches with the  position of the maximal negative shear stress in Fig. 10.
The erosion and deposition of sand leads to a temporal change in the height of the dune surface. In the steady state the dune moves shape invariantly and consequently the relation between changing sand flow rate and the deformation of the dune shape becomes constant. From the conservation of mass a term for the migration velocity of dunes (v com in our experiment) can be derived as [4,56] v com ∼ dI x,total dh = const.
In field measurements this simple formula is used to check the present steady state of a desert dune [56]. For our model dune the relation of sand flow rate I x,total (x) and dune height h(x) at a certain horizontal position x is plotted in Fig. 13. The dune is partitioned at the crest in an upstream side and a downstream side. Above a height of about 4 mm (two grain diameters) the slope dI x,total /dh becomes constant for both sides, an additional confirmation for the steady state. As mentioned above the bulk of the dune is composed of resting particles and creeping particles. To separate them, we use the threshold velocity v th = 5 mm/s as criteria. The value of v th is given by the standard deviation of the velocity distribution of all beads. Note, v th ≈ 1/10v com . It can be seen in Fig. 14 that sometimes beads in the core of the dune are set in motion, which results from the regrouping of the beads due to the erosion and deposition processes. Moreover, the snapshot in Fig. 14 illustrates that only a few particles are completely detached from the surface. This corresponds to ϕ v (x, y) shown in Fig. 8.
These observations are quantitatively shown in Fig. 15(a) with the distributions of the particle densities n r (resting), n c (creeping), and n s (saltating). For the analysis every recorded image is divided in columns with a width of 1 cm. The particle densities are determined as the time averaged number of particles counted in each column between t = 4 s and t = 5 s. For technical reasons the frames are shifted according to the center of mass of the dune.
The distributions of n r and n c in Fig. 15(a) show an excess of creeping particles on the downstream side, where the deposition occurs. At the upstream bottom of the dune n c and n s increase. They take constant values in the middle section and decrease in the deposition area. Note that the local densities n c and n s show that there are more creeping particles than saltating particles, but Fig. 12(b) shows that the flow rate of the saltating beads contribute the larger part of I x,total (x). This can be explained by the higher speed of the saltating particles.
The particle tracking method (see Fig. 4) allows for the investigation of transitions between the three types of beads in our experiment. Because every transition has two directions, we have to distinguish between rest-creep, creep-rest, creep-saltation, saltation-creep, rest-saltation, and saltation-rest. The number of transitions per second gives the transition rate. To simplify the notation, we label each connected pair with one variable: T rc (rest and creep), T cs (creep and saltation), and T rs (rest and saltation). The local distributions of the transition rates corresponding to Fig. 15(a) are shown in Fig. 15(b)-(d). In the areas of erosion and deposition T rc takes a local maximum. The maximum of T cs lies at the crest. The high values of all transition rates indicate the frequent collisions between the beads during the dune migration. The distribution of T rs in Fig. 15(d) shows that no reptation in the aeolian sense occurs in our experiment. The saltating beads do not originate from resting beads. They originate from creeping beads dragged from the fluid. The peak at the downstream side results from single beads, which are temporally not connected to the bulk of the dune. The net rate of erosion and deposition is obtained by subtracting the two parts of T rc , which yields ∆T rc in Fig. 15(e). The resulting curve fits to the erosion rate shown in Fig. 12(c). The distribution of n r in panel (a) reflects the contour of the bulk, shown in Fig. 15(f).
To conclude, we designed an experiment which allows the direct measurement of the particle dynamics and the surrounding driving water flow above and inside a twodimensional barchan dune model. The measurement of the particle dynamics inside the dune reveals the nature of the crawling motion of dunes, whose speed is one order of magnitude smaller than that of the creep and saltation of individual grains. This ratio is of course not universal, but specific for our downsized model dune. Due to the possibility to distinguish between creeping and saltating particles, we show that the erosion rate consists of comparable contributions from both. The saltation flow rate is slightly larger, whereas the number of saltating particles is one order of magnitude lower than that of the creeping ones. The velocity field of the saltating particles is comparable to the velocity field of the driving fluid, although that of the particles lags behind.
Moreover, it is demonstrated that the initial triangular heap evolves to a steady state with constant mass, shape, velocity, and packing fraction after one turnover time has elapsed. Within that time the mean distance between particles initially in contact reaches a value of approximately one quarter of the dune basis length, a number which serves to quantify the mixing process.
The shear stress of the water flow is shown to be related to the particle erosion rate, as well as the slope of the dune surface. It can be observed that the spatial profile of the shear stress reaches its maximum value upstream of the crest, while its minimum lies at the downstream foot of the dune. The particle tracking method reveals that the deposition of entrained particles occurs primarily in the region between these two extrema of the shear stress. Due to this deposition most of the mass of the dune is preserved. The upstream shift of the shear stress has been a crucial element for explaining the steady state in the minimal models. For a more fundamental understanding of the dune dynamics the presented facts are believed to encourage particulate models and simulations.