A Computer Simulation of Sand
Ripple Formation
By Peter Lamb, Jordan Kwan and Sam Ahn
May 1, 2002
1
Abstract
In recent years, scientists have begun to model sand dunes and ripples using
computer simulations. Building upon several papers we
constructed a computer model
which simulates the formation of sand ripples. We examined how the features of the
ripples were dependent on a variety of parameters. Our model was able to successfully
model the formation of primary and secondary ripples wi
th coarsening. We discovered
that the wavelength of the ripples were dependent on the wind velocity, but not creep due
to gravity or angle of repose.
Introduction
In desert regions around the world, one can commonly see beautiful sand patterns
which d
isplay remarkable similarity from desert to desert. Sand ripples are small, wave

like patterns which form in the sand and can range in wavelength from centimeters to feet.
Sand dunes are larger

scale piles of sand, in a variety of shapes. These patterns
have been
studied extensively in the twentieth century, intriguing researchers as to how and why they
form.
Two main effects, saltation and surface creep, are known to cause the formation of
the sand ripples and dunes. Saltation is the blowing of sand
grains from one location to
another by the wind. Surface creep is the shifting of the sand at the surface due to both the
wind and gravity. With the aid of computers, researchers are now able to better understand
sand dune and ripple formation. Computer
simulations can now closely model the creation
of these sand formations, and can be used to answer questions previously left unanswered.
One of the earlier attempts to simulate the formation of sand ripples was by
Nishimori and Ouchi. Using a very simp
lified cellular automata model, they implemented
the processes of saltation and creep due to gravity. They were able to show the virtual
formation of sand ripples, whose resulting wavelength was linearly proportional to the
strength of the wind, which is
in agreement with experimentation and observation of actual
sand ripples. They also presented a phase diagram showing whether ripples form under
varying parameters of wind and gravity.
A recent paper by Miao, Mu and Wu takes a cellular automata approach s
imilar to
Nishimori and Ouchi’s, but greatly expands on the model. In their model they include a
variety of other effects including avalanching on steep slopes, eddying on the downwind
sides of slopes, and saltation dependence on the local slope. Like Ni
shimori and Ouchi
they demonstrate the formation of ripples as well as the formation of a larger

scale sand
dune. They discuss the dependence of the ripple wavelength on the average hop distance.
In this project we successfully implement the simulation pr
esented in the paper by
Miao et al, with a few variations. We go beyond that which is presented in the paper by
examining the dependence of ripple formation on a variety of parameters.
2
Modeling Sand Ripples
As mentioned above, the majority of this mo
del is taken from those presented in
Nishimori, et al. and Miao, et al. The purpose of this project was to implement this model
on our own, examine the ripples’ dependence upon a variety of parameters not discussed in
either paper, and add what we could t
o the model. The equations we implemented are
presented below, most of which come from Nishimori and Ouchi. We will try to point out
which equations are different, and why we made the changes and additions we did. All of
the calculations detailed below
were applied in a simultaneous time step to every entry or
“square” in a 100 by 100 matrix. Each entry in the matrix had a numerical value which
represented the height of the sand in that square. A variety of parameters were entered into
the program desc
ribing wind strength, gravity, grain size, initial conditions, etc. To
conserve grains, those grains which blew off one edge entered on the other, and the edges
“communicate” with each other. Initial conditions varied, but for most of the calculations
a
normal random distribution centered at a height of seven or ten was used. The exact code
which generated our results was written in Matlab, and can be found in Appendix A.
There are three main elementary processes that occur in the desert due to wind.
They are
suspension
,
saltation
, and
sand creep
.
Suspension is when the sand grains are carried up in the air by the wind and cease
to fall back to the ground. One can imagine it like a bunch of particles in some liquid
medium, being swirled around and
suspended
in the medium. In this sense, sands that are
too light to fall back down are the particles and the air is the medium. Since suspension
keeps sand grains in the air, it is not involved in the formation of the patterns, and so is not
included in
either of the published models or our own modified version presented below.
Saltation happens when the sand grains are picked up by the wind and deposited
some distance later. Because the surface is not completely smooth, the wind will pick up
some of t
he grains that “stick out” and carry them some distance that is proportional to the
wind velocity. It is perhaps akin to sprinkling
salt
on some surface or substance.
The equations describing saltation are the following:
)
,
(
)
,
(
)
,
(
'
y
x
q
y
x
h
y
x
h
s
(1)
)
,
(
)
,
(
)
,
(
'
y
x
q
l
y
l
x
h
l
y
l
x
h
s
sy
sx
sy
sx
(2)
j
y
x
l
i
y
x
l
y
x
l
sy
sx
s
)
,
(
)
,
(
)
,
(
(3)
)))
,
(
tanh(
1
(
)
,
(
0
y
x
g
q
y
x
q
s
(4)
)
,
(
y
x
h
and
)
,
(
'
y
x
h
= heights of sand surface at point (x,y) before and after saltation
movement, respectively.
)
,
(
y
x
q
s
= transferred heights moving one step.
)
,
(
y
x
l
s
= horizontal displacement vector.
3
q0 is the average sand grain size,
)
,
(
y
x
q
s
is a function of the gradient of sand surface at
point (x,y), defined as:
2
2
)
,
(
y
h
x
h
x
h
sign
y
x
g
(5)
The components in (3) are defined as
))
tanh(
1
)(
(
)
,
(
0
x
h
h
w
l
y
x
l
x
x
sx
(6)
))
tanh(
1
)(
(
)
,
(
0
y
h
h
w
l
y
x
l
y
y
sy
(7)
where
x
l
0
and
y
l
0
represent the base distance hopped by all grains in the x and y directions
which is correlated with wind strength,
x
w
and
y
w
are other wind strength factors, and
x
h
and
y
h
are the slopes in x and y directions.
Equation (3) describes how far the sand moves in the x and y direction. Equations
(6) and (7) say that the saltation distance, the distance a sand grain “hops”, is directly
related to the wind strength and height of
the grain, and inversely related to the local slope.
In other words, the harder the wind is blowing, the farther the sand will fly. Also, the
higher a grain lies, the longer “airtime” it will have when it hops, and so will fly farther.
The
))
tanh(
1
(
y
h
factor says that on the upward wind

facing side of a developing ripple,
a saltating grain will fly less far because it is more likely to run into the face as it’s
jumping. Conversely, on the sheltered, lee

face, the grains have lots of open air
in front of
them and will fly farther.
Equation (4) describes how much sand gets lifted up at a time. There is a base
amount, q0, which is factored by
)))
,
(
tanh(
1
(
y
x
g
, which causes more sand to be
transported if the sand is on the windward slope,
and vice versa if on the leeward slope.
Sand creep is the general notion of sand being shifted along, or
creeping
on the
sand surface. There are two kinds: creep due to saltation and creep due to gravity impacts.
Creep due to saltation describes the
forward shifting of sand resulting from saltation
impacts. Compared to the other processes, it does not have a significant effect on the
ripple formation and so is not included in our model. Creep due to gravity is the settling
and flattening out of the
sand features by gravity, and can be described by:
NNN
n
n
NN
n
n
n
y
x
h
y
x
h
y
x
h
D
y
x
h
y
x
h
)
,
(
)
,
(
12
1
)
,
(
6
1
)
,
(
)
,
(
'
'
'
'
1
(8)
4
where
NN
n
y
x
h
)
,
(
'
represents the summation of the heights of the adjacently horizontal and
vertical sites and
NNN
n
y
x
h
)
,
(
'
represents the summation over
the sites diagonal to the site in
question. D is the rate of relaxation, which can be viewed as the strength of gravity. What
(8) says is that the height of a square in the next time step will differ by some proportional
amount D depending on the differ
ence between the weighted

average of the sites
surrounding the site and the site itself. Thus, if the site in question is overall higher than
the neighboring sites, then that height is decreased by some proportional amount, and vice
versa.
At some criti
cal angle, any sand which is piling up in a slope will “avalanche”
down until the critical angle is no longer exceeded. To model this phenomenon, we use
the following equations:
Applied to the lee slope:
x
x
y
x
h
y
y
x
h
y
)
,
(
)
,
(
(9)
if
tan
)
/(
))
,
(
)
,
(
(
2
/
1
2
2
y
x
y
y
x
x
h
y
x
h
, then let
2
/
)
tan
)
(
)
,
(
)
,
(
(
2
/
1
2
2
y
x
y
y
x
x
h
y
x
h
(10)
)
,
(
)
,
(
'
y
x
h
y
x
h
(11)
)
,
(
)
,
(
'
y
y
x
x
h
y
y
x
x
h
(12)
What this basically says is that if the slope in any direction of the site in question
exceeds a critical angle
, then in accordance to the steepness of the slope, a correction
factor
is measured that “avalanches” the site so that the resulting slope equals
.
Modeling Challenges
A simple model was first
constructed to see how accurate we could be
implementing just a few basic rules. To that end, we first implemented only the saltation.
What we found was complete randomness that looked like noise rather than ripples.
Adding the sand creep due to gravity
, however, corrected much of this. Since these
calculations take in to consideration the neighboring sites, in essence this process
connected all the sites together, and effectively smoothed out the randomness. With just
these two we were able to achieve
surfaces that looked like sand ripples. We further
refined our model, as described above, by adding in a gradient factor term to the saltation
lengths and grain transport amounts.
However, over time the ripples turned into spikes that realistically wou
ld not
happen. Therefore we attempted to implement the slope checking algorithm manifested in
equations (9) through (12) as described in Nishimori et al. Implementing this aspect of the
model was probably the most challenging aspect of the project.
One
large problem arose when we could not decide what should happen if there
was more than one direction whose slope exceeded the critical angle. Both directions had
5
to have slopes that fell to the critical angle, and that involved moving the sand to
appropr
iate amounts. The question became then, how much sand gets moved to each
square? Intuitively the answer should be some amount proportional to each slope, but
writing the code for this depends on how many neighboring sites actually need to be
adjusted, wh
ich quickly becomes very complicated.
Another problem had to do with over counting. For each site, there are 8
neighboring sites, the four “immediate” neighbors to the sides and the four “secondary”
diagonal neighbors. For each time step, each site nee
ds to be checked with all of its
neighboring sites once and only once. If one site is checked against all eight of its
neighboring sites, then double counting occurs when the neighbor himself checks. We had
to find a way to correct for this double counti
ng problem along with the sand movement
problem.
To solve the latter problem, we decided to check only four of the eight neighboring
sites, per every site in question. We checked all three squares below and the square to the
immediate right. If applied
to every square, then all directions will eventually be checked,
without over

counting. Then, we identified which directions had overly steep slopes, and
created a system of equations that, when solved, determined how much sand needed to be
adjusted to e
ach square.
Verification of Numerical Techniques
To verify the validity of our model, we conducted several tests. When run with no
wind, the result is that the grid gets flattened out. No ripples form, and the irregularities in
the initial grid get sm
oothed due to creep from gravity. Eliminating the creep from gravity
would result in no movement of the sand. This makes physical sense.
Our model also conserves mass. That is, any grain that is removed from a
particular cell gets added to another cell.
Over time, the sum of the grains over all cells
remains constant. In the real world, grains of sand do not disappear over time when wind
blows over them. They are just moved. Therefore, we are pleased with the conservation
of mass in our model. Energ
y is not conserved because energy is constantly being added
into the system from the wind.
Finally, our model is stable provided that our critical angle is reasonable and there
is reasonable creep from gravity. Given these conditions, our resulting grid d
oes not
display any awkward spikes that would not be physically possible. However, without an
angle of repose and without sufficient creep due to gravity, our model displays an
instability where the ripples keep the same wavelength and with the amplitude
increasing
to unreasonable heights. In these cases, it is clear that our model falls apart, though we are
not concerned with this since there is no analogous case in the real world.
A concern that we do have is that our model wraps the grid around. When
a grain
of sand flies off one end of the grid, it blows in from the other end. We believe that this
might affect the wavelength of the ripples. Furthermore, we are concerned that this may
drastically affect our results when the grains blow further than t
he length of our grid
(blows further than n cells when our grid has length c). Thus, we tried to avoid these cases
by keeping the wind sufficiently low. Naturally, grains of sand do not wrap around in the
real world. However, we were forced to implement
our model in this way to simulate a
large field of sand since computationally, it was only feasible to model a relatively small
6
100 by 100 grid. If we were able to increase the grid size sufficiently, these wraparound
effects would be reduced, and we mig
ht be able to eliminate the wraparound entirely.
Another concern is that our model is not entirely symmetric. In nature, the
direction in which the wind is blowing does not matter. The ripples will always align
themselves perpendicular to the direction
of the wind. Initially, our model could only have
wind blowing in the x

direction. Later we added a y component to the wind, so that ripples
formed very nicely in either the y or x directions. However, because a matrix does not
look “square” when viewed
at a diagonal, as might be expected the ripples have a hard
time settling in and “lining up” when the wind blows at a diagonal. They do eventually
form, but they are much less stable than those in the purely x or y directions and take much
longer to form
.
It should be noted that since our model did not take into account any sort of units,
these results are primarily qualitative. The model is useful for measuring relative effects.
Results
After getting the model working well, we began testing the sand
ripples
dependence on a variety of parameters. Figure 1 shows some images from a typical sand
ripple simulation where the ripples coarsen into larger ripples.
(a)
(b)
7
(c)
(d)
Figure 1

Sand ripple simulation with
hopX=30, windX=0.5, grain=0.1, gravity=0.9
shown after (a) 30 cycles (b) 60 cycles (c) 120 cycles and (d) 180 cycles. This is on a 100
by 100 grid, with initial conditions distributed normally with mean of 10. The initial
wavelength is about 10 (hopX)
, before coarsening into a secondary wavelength of around
25.
Observation and experimental results indicate that the wavelength of the ripples is
dependent upon the average saltation length (Bagnold and Lancaster). This depends
primarily upon the streng
th of the wind. The other sand ripple models show this
dependence. Using our model we also decided to see how sand ripple wavelength depends
upon our wind strength parameter and our mean hop length parameter.
Figure 2 shows how the primary wavelength
and the secondary wavelength after
the ripples coarsen increases as we increase windX. The primary wavelength increases
proportionally as expected. There are a few points which decrease with increasing
wavelength. We’re not quite sure why this is, but o
bviously it is some defect in our model.
More interesting is what is occurring with the secondary wavelength. The coarsening is
nonexistent, then secondary ripples appear, then disappear, then reappear again. We’re not
sure why this happens, but is some
thing to explore in the future.
8
0
1
2
3
4
0
50
100
Ripple Wavelengt h vs. Wind Velocity
wavelength
subwavelength
windX
Figure 2

Wavelength dependence upon windX, our wind speed parameter, after 100 time
cycles. HopX = 20, grain size = .1, gravity = .95
Figure 3 shows how wavelength scales with our base grain hop,
which is what was
explored in the other papers. This more clearly has a direct correlation, with less weird
points.
0
10
20
30
40
0
50
100
Ripple Wavelengt h vs. Base Grain Hop
55
0
Wavelength
SubWavelength
38
0
HopX
Figure 3

Wavelength dependence upon hopX after 100 time cycles. WindX=0.5, grain
size = 0.1, gravity = 0.95
Tests
were also run that examined the effects of the gravity parameter. Figure 4
shows an example of the results of varying the creep due to gravity. Since the creep due to
gravity smoothes out large spikes in the sand, having a low gravity parameter resulted
in
graphs that tended to be very spiky. In fact, setting that parameter to zero resulted in what
9
appeared to be random noise. As you can see in figure 4, increasing the gravity parameter
made the resulting graph smoother. We further discovered through
testing that adjusting
the gravity parameter had no effect on the wavelength (as this parameter did not affect the
hop distance of the grains). Surprisingly, the gravity parameter had no significant effect on
the amplitude of the ripples aside from leveli
ng the most extreme spikes.
(a)
(b)
(c)
Figure 4
–
Ripple dependence on gravity for hopX=10, windX=0.5, grain = 0.1 for 60
time steps. (a) gravity = 0.3 (b) gravity = 0.7 (c) gravity = 1.0
10
Our program also accoun
ted for how the sand cannot achieve a slope higher than
some critical angle, due to physically intuitive reasons. The effects of this self

reposing
mechanism can be seen in Figures 5 and 6. Figures 5 shows what the ripples might look
like without the sel
f

reposing mechanism. Contrasting this with Figure 6, one can see that
the ripples form more skewed shapes that more closely reflect sand ripples with the wind
blowing on one side. The slopes on the leeward side fall sharply away from its peak,
which mir
rors how the wind pushes against the ripples on one side.
Figure 5:
A cross section of results, showing a profile of ripples without the self

reposing
mechanism
Figure 6:
The same simulation run with exact parameters as in Figure 5, except now
con
taining the self

reposing mechanism
Conclusions
Figure 1 illustrates the typical evolution of the computer simulation. This reflects
how the computer simulation process approaches a steady state, creating consistent,
relatively uniform ripples that run
perpendicular to the direction of wind, regardless of the
initial conditions. Figures 2 and 3 show a clear relationship between the saltation
parameters and ripple wavelength. Consistent with research and empirical data, ripple
wavelengths in our model
are proportional to the wind velocity. Figure 4 shows the effects
of gravity. As was expected, larger gravity parameters smooth out the surface more.
Qualitatively, the amplitude of the waves is only affected by gravity in so much as each
ripple is smoo
thed out. Consequently, some sharp ripples have their heights shortened
because of the smoothing out. Additionally, the self

reposing mechanism was successful,
in that it models how the ripples “lean” more in the direction of the wind. Interestingly,
ho
wever, the self

reposing mechanism does not seem to affect the wavelength or the
amplitude. Overall, akin to existing literature, the formation of sand ripples has been
shown to be a self

organizing process.
11
References
Nishimori and Ouchi, “Formation
of Ripple Patterns and Dunes by Wind

Blown Sand,”
Physical Review Letters
, The American Physical Society, 1993.
Miao, Mu, and Wu, “Computer Simulation of Aeolian Sand Ripples and Dunes,”
Physics
Letters A,
Elsevier Scince B.V., 2001.
Bagnold, R. A.
The Physics of Blown Sand and Desert Dunes
. William Morrow & Co.:
New York, 1941.
Lancaster, Nicholas.
Geomorphology of Desert Dunes.
Routledge: New York, 1995
Appendix A
–
Sand ripple model MatLab m

file code
function [Hnew,M] = ripples2(H,hopX,wi
ndX,hopY,windY,grain,gravity,critAng,numsteps)
[rows cols] = size(H);
Heven = H;
Hodd = H;
for currstep = 1:numsteps
% Hodd blowing to Heven
Heven = Hodd;
for x = 1:rows
for y = 1:cols
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% here we define the coordinates for all surrounding neighbors of a
% point
if (x==1)
xUp=rows;
else
xUp = x

1;
end
if (x==rows)
xDown=1;
else
xDown = x+1;
end
if (y==1)
yLeft=cols;
else
yLeft = y

1;
end
12
if (y==cols)
yRt=1;
else
yRt = y+1;
end
h = Hodd(x,y);
hR = Hodd(x,yRt); hRD = Hodd(xDown,yRt); hRU = Hodd(xUp,yRt);
hL = Hodd(x,yLeft); hLD = Hodd(xDown,yLeft); hLU = Hodd(xUp,yLeft);
hD = Hodd(xDown,y); hU = Hodd(xUp,y);
%%%%%%%%%%%%%%%%%%
% SALTATION
% we define the gradients
delHx = Hodd(xDown,y)

Hodd(x,y);
delHy = Hodd(x,yRt)

Hodd(x,y);
delH = sign(delHx)*sqrt(delHx^2 + delHy^2);
% the hop length depends on the height and slope
hopLengthX = (hopX +
windX*Hodd(x,y))*(1

tanh(delHx));
hopLengthY = (hopY + windY*Hodd(x,y))*(1

tanh(delHy));
% the amount of grains transported depends on the slope, delH
grainAmt = (grain)*(1+tanh(delH));
% this is where the grains blow
blowToX = round(x + hopLengthX);
blowToY = round(y + hopLengthY);
if (blowToX > rows)
blowToX = blowToX

rows;
en
d
if (blowToY > cols)
blowToY = blowToY

rows;
end
Heven(x,y) = Heven(x,y)

grainAmt;
Heven(blowToX,blowToY) = Heven(blowToX,blowToY) + grainAmt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CREEP due to gravity
% here we're taking a weighted sum of the 4 immediate neighbors
FirstNbrSum = 1/6*( hU + hD + hL + hR );
% we take a less

weighted sum of the four "corner neighbors"
SecondNbrSum = 1/12*( hLU + hRU + hLD + hRD);
Heven(x,y) = Heven(x,y) + gravity*(FirstNbrSum+SecondNbrSum

h);
13
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AVALANCHE if the slope gets too big
sandIndex = [];
b = [];
if ((h

hR)/1 > tan(critAng))
sandIndex = [sa
ndIndex;1];
b = [b;tan(critAng)+hR

h];
end
if (h

hRD)/(sqrt(2)) > tan(critAng)
sandIndex = [sandIndex;2];
b = [b;sqrt(2)*tan(critAng)+hRD

h];
end
if (h

hD)/1 > tan(critAng)
sandIndex = [sandIndex;3];
b = [b;tan(critAng)+hD

h];
end
if (h

hLD)/(sqrt(2)) > tan(critAng)
sandIndex = [sandIndex;4];
b = [b;sqrt(2)*tan(critAng)+hLD

h];
end
n = length(b);
A =

ones(n);
for i = 1:n
A(i,i)=

2;
end
sandShift = A
\
b;
for j = 1:n
if sandIndex(j) == 1
Heven(x,yRt)=hR + sandShift(j);
elseif sandIndex(j) == 2
Heven(xDown,yRt)=hRD + sandShift(j);
elseif sandIndex(j) == 3
Heven(xDown,y)=hD + sandShift(j);
elseif sandIndex(j) == 4
Heven(xDown,yLeft)=hLD + sandShift(j);
end
end
Heven(x,y)

sum(sandShift);
end
end
Hodd = Heven;
meshz(Hodd)
axis([0 100 0 100 0 50])
pause(.01)
14
if mod(currstep,3)==0
M(currstep/3)=getframe;
end
end
Hnew = Hodd;
Comments 0
Log in to post a comment