# A Particle Model of Rolling Grain Ripples Under Waves

## I The rolling grain ripples

An example of a flat bed with rolling grain ripples coexisting with vortex ripples is seen in Fig. 1. In the middle of the picture the ripples have not yet formed and the bed is still flat, while on the top vortex ripples are seen to invade into the flat bed. The vortex ripples are typically nucleated from the boundaries or from a pertubation in the bed. In the lower part of the picture the rolling grains ripples have formed on the flat bed, and are seen as the small bands of loose grains on top of the flat bed.

The flow over the bed created by the surface wave is oscillating back and forth in a harmonic fashion. This flow creates a shear stress on the bed $`\tau (t)`$ which in non-dimensional form reads:

$$\theta (t)=\frac{\tau (t)}{\rho (s1)gd},$$
(1)   
where $`\rho `$ is the density of water, $`s`$ is the relative density of the sand (for quartz sand in water $`s=2.65`$), $`g`$ is the gravity and $`d`$ is the mean diameter of the grains. $`\theta `$ is usually called the Shields parameter . When the shear stress exceeds a critical value $`\theta _c`$, the grains start to move. For a turbulent boundary layer the value of the critical Shields parameter is $`\theta _c0.06`$ . The grains which have become loosened from the bed, start to move back and forth on the flat bed, and after a while the grains come to rest in parallel bands. In the lee side of each band the bed is shielded from the full force of the flow, creating a “shadow zone” where the grains move slower than in the upstream side of the bands.

Due to this shadow zone, more grains end up in the bands than leave the bands, and they grow until they form small ridges, the rolling grain ripples. When the rolling grain ripples are fully developed, no grains will be pulled loose from the bed in the space between them, and they are stable. However, in reality the rolling grain ripples are dominated by invading vortex ripples, which is the main reason why they are rarely observed in Nature.

## II A simple model

The above scenario can be formulated mathematically by writing an equation of motion for each grain/particle. In the following an equation of motion for the particles are developed. In the beginning the particles represents the grains, but as the single grains quickly merge, the particles most of the time represents a ripple. First the velocity of each particle is found, assuming that the particle is alone of the flat bed, and then the influence of the shadow zones from neighboring particles are taken into account.

Consider $`N`$ particles rolling on top of a rough, solid surface. Each particle is characterized by its position $`x_i`$ and its height $`h_i`$ (see Fig. 2). As the ripples are triangular, the area of each particle $`A_i`$ and their heights are related as:

$$h_i=\sqrt{A_i\mathrm{tan}\varphi },$$
(2)   
where $`\varphi `$ is the angle of repose of the sand (33).

The ripple moves back and forth slower than a single grain, according to the “1/height” law. This law is well known in the study of dunes in the desert  or sub-aqueous dunes , and can be illustrated by a simple geometrical argument. Suppose there is a flux of sand over the crest of a ripple or a dune $`q_{crest}`$ (Fig. 3). To make the ripple move a distance $`\delta x`$ an amount of sand $`h\delta x`$ is needed. As the sediment flux is amount of sand per unit time, the velocity of the ripple is $`u_{ripple}=q_{crest}/h1/h`$. If the height of the initial particles (the single grains) is assumed to be equal to the grain diameter $`d`$, the velocity of the particles can be related to that of the single grains as:

$$u_i=\frac{d}{h_i}U_g\mathrm{sin}(\omega t),$$
(3)   
where $`U_g`$ is the velocity amplitude of the motion of a single grain and $`\omega `$ is the angular frequency of the oscillatory motion.

The flow moves back and forth and makes the particles move accordingly on top of the bed. In the wake of each particle/ripple there is a shadow zone (Fig. 4), which is the area behind the particle where the absolute value of the shear stress is smaller than it would be on a flat bed. The length of the shadow zone is therefore larger than the length of the separation bubble formed by the particle (note that the shadow zone would be present even in the absence of separation). If the shadow zone is much smaller than the amplitude of water motion, the flow in the lee side of the ripple can be assumed to have sufficient time to become fully developed. The fully developed flow over a triangle is similar to that past a backward facing step in steady flow, which has been extensively studied (see e.g. Tjerry ). In that case the relevant quantities, i.e., the length of the separation bubble, the length of the wake etc., scale with the height of the step. As a first assumption, the shadow zone is therefore assumed to have a length which is proportional to the height of the particle: $`\alpha _sh_i`$. If a particle enters the shadow zone of another particle, it is slowed down according to the distance between the particles. This means that the actual velocity of a particle is $`u_if(\mathrm{\Delta }x)`$ where $`\mathrm{\Delta }x`$ is the distance between the grain and the nearest neighbor upstream, $`u_i`$ is the velocity of the particle outside the shadow zone and $`f`$ is a function determining the nature of the slowing of the motion of the particle. A simple linear function is used, as shown in Fig. 4. The exact form of the function $`f`$ is not crucial, as will be evident later, the important parameter being the extent of the shadow zone as determined by $`\alpha _s`$.

It is now possible to write the equations of motion for the particles as a system of coupled ODEs:

$$\dot{x_i}=\frac{d}{h_i}u_g(t)\underset{\mathrm{positive}\mathrm{half}\mathrm{period}}{\underset{}{f\left(\frac{x_ix_{i1}}{\alpha _sh_{i1}}\frac{u(t)}{|u(t)|}\right)}}\underset{\mathrm{negative}\mathrm{half}\mathrm{period}}{\underset{}{f\left(\frac{x_ix_{i+1}}{\alpha _sh_{i+1}}\frac{u(t)}{|u(t)|}\right)}}$$
(4)   
for $`i=1,\mathrm{},N`$, where $`u_g(t)=U_g\mathrm{sin}(\omega t)`$. The motion of a particle is thus made up of three parts (i) the motion of the single undisturbed particle, (ii) the effect of the shadow from the particle to the left $`(i1)`$, which might affect particle $`i`$ in the first half period and (iii) the effect of the shadow of the particle to the right $`(i+1)`$ in the second half period.

When lengths are scaled by the diameter of the grains and time by the frequency $`\omega `$, it is possible to identify the three relevant dimensionless parameters of the model:

$`\alpha _s`$, the length of the shadow zone divided by the height.

$`a_g/d`$, the amplitude of the motion of a single grain, divided by the grain diameter ($`a_g=\omega U_g`$).

$`\lambda _i/d`$, the initial distance between the grains; $`\lambda _i=L/N`$, where $`L`$ is the length of the domain and $`N`$ is the initial number of grains.

### A Relation to physical quantities

Even though the model seem quite heuristic, the parameters entering the model, $`a_g/d`$, $`\lambda _i/d`$ and $`\alpha _s`$ can be related to physical parameters describing the flow and the properties of the grains. The line of arguments presented here closely follows those used to derive the flux of sand on a flat bed (the bed load), as can be found in text books, i.e., , or in . First the velocity of a single grain will be derived, from which $`a_g/d`$ can be inferred. Thereafter the initial number of grains in motion is found, from which follows $`\lambda _i/d`$. Finally the length of the shadow zone $`\alpha _s`$ is discussed.

#### 1 The velocity of the grains

The velocity of the grain can be found by considering the force balance on a single grain lying on the flat bed. The grain is subject to a drag force proportional to the square of the relative flow velocity $`u_r=u_{nb}u_g`$ where $`u_{nb}`$ is the velocity near the bed and $`u_g`$ is the velocity of the grain:

$$F_d=\frac{1}{2}C_D\rho A|u_r|u_r,$$
(5)   
where $`A`$ is the area of the grain and $`C_D`$ is a drag coefficient. The numerical signs is used to obtain the right sign of the force. The velocity profile in the vicinity of the bed is supposed to be logarithmic. As shown in  this is a reasonable assumption except when the flow reverses. However, during reversal the velocities are anyway small and the accuracy is of minor importance. The logarithmic profile over a rough bed can be written as (eg. ):

$$u(y)=\frac{u_f}{\kappa }\mathrm{ln}\left(\frac{30y}{k_N}\right),$$
(6)   
where $`\kappa =0.41`$ is the von Kármán constant, $`u_f\sqrt{\tau /\rho }`$ is the friction velocity and $`k_N`$ is the Nikuradse roughness length. It is then possible to find the near bed velocity as the velocity at $`y1/2d`$:

$$u_{nb}=\xi u_f$$
(7)   
where the constant $`\xi `$ can be determined from Eq. (6) by assuming $`k_N=d`$. Opposing the drag on the grain is the friction of the bed:

$$F_f=\mu W$$
(8)   
where $`\mu `$ is a friction coefficient and $`W=\rho g(s1)d^3\pi /6`$ is the immersed weight of the grain. By making a balance of forces, $`F_d+F_f=0`$, the velocity of the grain can be found:

$$u_g=\xi u_f\left(1\sqrt{\left|\frac{\theta _c}{\theta }\right|}\right),$$
(9)   
where

$$\theta _c=\frac{4\mu }{3C_D\xi ^2},$$
(10)   
is the critical Shields parameter. Usually $`\mu =\mathrm{tan}\varphi 0.65`$ .

#### 2 The initial spacing of grains

Now that the velocity of the grains has been calculated, there still remains to determine the number of grains per area $`n`$ in motion. To this end a small volume of moving sand at the top the flat bed is considered. The balance of the forces acting on this volume is written as:

$$\tau _b=\tau _G+\tau _c.$$
(11)   
The interpretation of the terms is as follows: The parameter $`\tau _b`$ is the shear stress on the top of the bed load layer. It is assumed that this is equal to the shear stress on a fixed flat bed. $`\tau _G`$ is the stress arising from the inter-granular collisions, giving rise to “grain-stresses” , modeled as: $`\tau _G=n\mu W`$. It is assumed that the inter-granular stress absorbs all the stress, except the critical stress $`\tau _c`$; this is the so-called “Bagnold hypothesis” . Making Eq. (11) non-dimensional by dividing with $`\rho (s1)gd`$, the number of grains in motion is found as:

$$n=\frac{6}{\pi d^2\mu }(\theta \theta _c),$$
(12)   
If $`\theta <\theta _c`$, then no grains are in motion and $`n=0`$. $`n`$ can also be viewed as the initial density of grains, and by assuming a square packing of the grains the initial distance between the grains becomes:

$`{\displaystyle \frac{\lambda _i}{d}}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{n}d}}`$ (13)   
$`=`$ $`\sqrt{{\displaystyle \frac{\pi \mu }{6(\theta \theta _c)}}}.`$ (14)   
#### 3 The length of the shadow zone

The last parameter $`\alpha _s`$, which characterize the length of the shadow zone, is a bit harder to estimate accurately. By exploring the analogy with the backward facing step which was suggested in section II, we can get some bounds for $`\alpha _s`$. In the backward facing step there is a zone with flow separation which extends approximately 6 step heights from the step, see Fig. 5. After approximately 16 step heights there is a point where the shear stress has a small maximum. It follow that the length of the shadow zone should be longer than the separation zone, but shorter that the point of the maximum in the shear stress, i.e., $`6<\alpha _s<16`$.

### B Numerical and analytical solutions of the model

In the following section the behavior of the model is examined. To study the detailed behavior the set of coupled ordinary differential equations (4) is integrated numerically. It is seen that the model reaches a steady state and an analytical expression for the spacing between the ripples in the steady state is developed.

The numerical simulations in this section are based on a simple example with $`a_g/d=35`$ and $`\lambda _i/d=3.23`$ and $`\alpha _s`$ is set to 10. As initial condition all particles have an area of $`1.0\pm `$ 10 %, to add some perturbation. The initial number of particles $`N`$ in this example is 800.

In the first few periods a lot of grains is colliding and merging (Fig. 6). As the ripples are formed and grow bigger the evolution slows down, until a steady state is reached (Fig. 7). The spacing between the ripples in the steady state show some scatter around the average value, which is also observed in experiments. The variation in the average spacing of the steady state $`\lambda _{eq}`$ for realizations with different initial random seed, turned out to be on the order of $`1/N_{eq}`$ where $`N_{eq}`$ is the final number of ripples.

To study the behavior of the average spacing of the ripples $`\lambda _{eq}`$, a number of simulations were made where the parameters were varied one at a time. Each run was started from the initial disordered state.

Changing $`a_g/d`$ only results in a minor change in the spacing of the ripples (Fig. 8a). The final spacing between the ripples does depend on the length of the shadow zone $`\alpha _s`$; the longer the shadow zone, the longer the ripples (Fig. 8b). This can be used to estimate the average equilibrium spacing between the ripples. When the distance between two ripples is longer than the shadow zone of the ripples, they are no longer able to interact. This gives:

$$\lambda _{eq}>\alpha _sh_{eq},$$
(15)   
where subscript $`eq`$ denotes average value at equilibrium. However if two ripples have a spacing just barely shorter than Eq. (15), they will be able to interact and eventually they will merge. One can therefore expect to find spacings up to $`\lambda _{eq}<2\alpha _sh_{eq}`$. Assuming that the average length is in between the two bounds one get:

$$\lambda _{eq}=\gamma \alpha _sh_{eq},1<\gamma <2.$$
(16)   
where $`\gamma `$ can be found by comparing the results from the full simulations with Eq. (16). The height of the ripples at equilibrium can be found by splitting the initial number of particles evenly onto the equilibrium ripples. Then the average area at equilibrium is $`A_{eq}=\lambda _{eq}/\lambda _id^2`$ and from Eq. (2) follows that the height is $`h_{eq}/d=\sqrt{\mathrm{tan}\varphi \lambda _{eq}/\lambda _i}`$, which gives an average equilibrium spacing:

$`{\displaystyle \frac{\lambda _{eq}}{d}}`$ $`=`$ $`\gamma ^2{\displaystyle \frac{\alpha _s^2d\mathrm{tan}\varphi }{\lambda _i}}`$ (17)   
$`=`$ $`\alpha _s^2\gamma ^2\sqrt{{\displaystyle \frac{6\mathrm{tan}\varphi }{\pi }}}\sqrt{\theta \theta _c}.`$ (18)   
The equilibrium spacing is therefore found to be proportional to $`\sqrt{\theta \theta _c}`$ with the constant of proportionality being made up of $`\alpha _s`$, $`\gamma `$ and various geometrical factors. All the quantities related to the dynamical evolution of the ripples, i.e., the velocity of the ripples, the shape of the function $`f(\mathrm{\Delta }x)`$ etc. does not enter into the expression.

### C Comparison with experiments

The only parameter that has not been accurately determined is $`\alpha _s`$. The value of this parameter can be estimated by comparison with measurements.

In 1976, Sleath made a series of experiments, measuring spacing between rolling grain ripples . The ripples were formed on a flat tray oscillating in still water using sand of two different grain sizes: 0.4 mm and 1.14 mm. To compare with the experiments the value of $`\theta `$ need to be calculated. $`\theta `$ reflect the number of grains in motion, and it seems natural to use the maximum value during the wave period, $`\theta _{max}`$. To estimate $`\theta _{max}`$ the shear stress on the bed has to be estimated. The maximum shear stress on the bed, $`\tau _{max}`$ during a wave period can be found using the concept of a constant friction factor $`f_w`$ ():

$$\tau _{max}=\frac{1}{2}\rho f_wu_{max}^2,$$
(19)   
with $`u_{max}`$ being the maximum near-bed velocity. The friction factor can be estimated using the empirical relation :

$$f_w=0.04\left(\frac{a}{k_N}\right)^{0.25}$$
(20)   
where $`a=u_{max}\omega `$ and $`k_Nd`$.

The range of Shields parameters in the experiments was found in this way to range from the critical Shields parameter to $`\theta =0.42`$. For the high Shields parameters the rolling grain ripples were reported to be very unstable, and to quickly develop into vortex ripples. In these cases, the measured ripple spacing then reflects the spacing between the rolling grain ripples before they developed into vortex ripples .

In Fig. 9 the experimental results are compared with runs of the model using $`\alpha _s=15.0`$ (the reason for this particular value will be clear shortly) and $`N=10000`$. By fitting all the runs to Eq. (16) it was found that $`\gamma =1.40`$. The results using Eq. (18) and $`\gamma =1.40`$ are shown with a line.

First of all, it is seen that Eq. (18) predicts the results from the full model Eq.(4) well. The correspondence between the model and the experiments is reasonable, but there are some systematic discrepancies, which will be discussed.

There are some few points with small ripple spacing for which the model does not fit the measurements. These measurements has a Shields parameter very near the critical (i.e., just around the onset of grain motion), which implies some additional complications. The grains used in the experiment were not of a uniform size; rather they were part of a distribution of grain sizes, and the grain size reported is then the median of the distribution, $`d_{50}`$. The Shields parameter is calculated using the median of the distribution, but actually one could calculate a Shields parameter for different fractions of the distribution, thus creating a $`\theta _{10}`$, a $`\theta _{50}`$ etc.. When $`\theta _{50}`$ is smaller than then the critical Shields parameter, $`\theta _{10}`$ might still be higher than the critical Shields parameter. This implies that grains with a diameter smaller than $`d_{50}`$ will be in motion, while the larger grains will stay in the bed. As only $`d_{50}`$ is used in the calculation of the equilibrium ripple spacing, the distance between the grains $`\lambda _i`$ will be overestimated near the critical Shields parameter, where the effect of the poly-dispersity is expected to be strongest. An overestimation of $`\lambda _i`$ will lead to an under-prediction of the ripple length, which is exactly what is seen in Fig. 9.

There are also three points from the experiments taken at very large Shields parameters, which seems to be a bit outside the prediction of the model. As already mentioned, these points are probably doubtful because of the very fast growth of vortex ripples. It is therefore reasonable to assume that vortex ripples invaded the rolling grain ripples before these had time to reach their full length.

To find a reasonable value of $`\alpha _s`$ Eq. (18) was fitted to the experimental points. To avoid the points which might be of doubtful quality, as discussed above, only the points in the range $`0.075<\theta <0.3`$ were used. This gave the value of $`\alpha _s=15.0`$, in agreement with the qualitative arguments in section II A 3.

## III Discussion of the results

From the comparison of the model with measurements it seems as if the model confidently reproduces the experiments.

In the model the number of grains in motion is constant (even though the number of particles changes). In an experimental situation, however, new grains might be lifted from the bed and added to the initial number of grains in motion. As the part of the flat bed between the ripples is covered by the shadow zones of the particles, these stretches will be shielded from the full force of the flow, and only very slowly will new grains be loosened here. This small addition of new grains will result in a slow growth of the rolling grain ripples, such that eventually grow into vortex ripples. This slow growth is very well illustrated by recent measurements , but not covered by the present model.

## IV Conclusion

In conclusion, a model has been created which explain the creation and the equilibrium state of rolling grain ripples of the type described by Bagnold. The final distance between the ripples is proportional to $`\sqrt{\theta \theta _c}`$. The model has been compared with measurements with reasonable agreement.

###### Acknowledgements.

It is a pleasure to thank Tomas Bohr, Clive Ellegaard, Enrico Foti, Jørgen Fredsøe and Vachtang Putkaradze for useful discussions. I also wish to thank the anonymous referees for constructive criticism.