Labs Sediment Transport Mechanics


Sediment Transport Mechanics by Greg Tucker & Bob Anderson, Spring 2013

Modeling assignment # 1. Experiments with water volume balance in buckets and basins

Goals: You will gain experience with: setting up a code in IRF format; setting up simple arrays; stepping incrementally in time; and in plotting your results.
a) Write a MATLAB code that answers the following question: given a fixed input of water to the bucket, and a hole in the base of the bucket, how will the water level in the bucket evolve? Assume that the bucket starts with no water in it. Produce a plot of water level as a function of time.

b) Now allow the inputs to vary through time, in particular allow the input to oscillate sinusoidally. Specify a period and amplitude of the oscillation about the mean. Again, answer how the water level in the tank will vary through time. Produce a two-part plot in which you display both input and output of water on a top plot, and the water level on the lower plot.

Answers to Leaky Bucket assignment # 1a &1b.

1.Bob Anderson's demonstration code

Closed Lake Basin # 1c &1d.


Now let’s take this to a little bit more real setting, a closed basin lake in which there is input solely from runoff from the surrounding landscape, and outputs solely from evaporation from the lake surface. One of the most famous examples in the US is Mono Lake in eastern California at the foot of the Sierras.

We start with a generic lake basin, and move on to a more real basin after developing a simple model. Consider a simplistic closed lake basin, called Snowcone Basin, shown in the figure below.

The lake basin: The drainage basin has an area Ab. There are no outlets. The geometry of the basin is a cone, with the tip of the cone at an elevation of 1000 m, and a spread angle for the cone of a.

The climate: On the average, it precipitates in the basin, uniformly (by which I mean evenly over the entire area), at a rate P (in meters of water equivalent per year). The climate is such that it also evaporates off any standing water in the basin at a mean annual rate E (also in meters per year).

The runoff: The precipitation is imperfectly converted to runoff when the precipitation falls on any surface but the lake itself, some of it being lost to transpiration and evaporation from the surface of the basin and its vegetation before the water can make it to any lake in the basin. The conversion factor, called the runoff coefficient, is . Of course, any precipitation that falls on the lake is not affected by this factor. Given this information, let us now ask several questions of the lake in this basin. Note that in all cases, E>P. (You should ponder how you could possibly get a lake under such circumstances!)

a) Assume that the lake in the basin has achieved a steady state condition, meaning that it has found an equilibrium in which the lake level and the lake area are constant in time. Write out the equation for the water balance that this implies (rate of volume into the lake = rate of volume out of the lake), and solve it for the expected altitude of the lake surface, z.

b) Now let’s ponder a specific basin. Let the area of the basin be 1000 km2 and the spread angle of the basin be 4°, an angle typical of the slopes of alluvial fans. Let the precipitation, P, be 0.1 m/yr, the runoff coefficient be 0.5, and the evaporation, E, be 2.5 m/yr. Assuming that these conditions have prevailed for some long time, at what altitude, zl, do you expect to see the modern shoreline? Now plot the expected lake level as a function of a range of possible mean annual precipitation, from P=0 to P=0.8m/yr.

c) We now travel to the edges of the modern lake in Snowcone Basin, and discover that there are lake shores well above the level of the modern lake. We measure the altitude of the paleo-shoreline to be 1250 m. Quantitatively, what can you tell us about the climatic conditions that might have prevailed at the time of the formation of this shoreline? I would like to know by what percentage of the modern precipitation the past precipitation must have been enhanced to raise the lake to this level, if the evaporation is held constant, and by how what percentage the evaporation in the past must have been reduced relative to the present, if the precipitation is held constant. (Likely the answer is of course some combination of these effects...and please assume the runoff coefficient is also unchanging.)

That gives you some sense of what controls the steady state lake level. Another type of question we can ask of this simple system is one of time scales. Think about how fast the system should respond to changes in the climatic forcing. While the basin geometry remains constant over these time scales, as it is determined by long term tectonic and landscape evolution processes, the climate can change very rapidly, causing changes in the lake surface elevation, the hence lake area and the lake volume. Clearly, this is important to think about if the lake water is being put to some human use, or if we have to deal in some way with expanded lake surfaces. For example, I-80 runs right along the shoreline of the Great Salt Lake, at an elevation not more than a few meters above it!

To answer this last question, construct a lake model in which you track water balance in the lake, and predict a history of lake level. We ask what is the characteristic time scale for responding to a change in climate. Does it take a decade or a century to respond? Is there hysteresis, meaning does it respond with one timescale when rising, another when lowering?

d) To accomplish this, construct a generic lake model, and spin it up from zero water in the basin until it reaches steady state (or very close to it). You will have to explore to find the appropriate time step and the appropriate model duration. Please show us histories of lake level, area and volume. Now assess this characteristic time of the lake. In radioactive systems we often use the “half-life” of a system, but more formally, it is better to define the characteristic time by determining the time it takes to reach (1-(1/e)) of its steady state value.

Answers to Closed Lake assignment # 1c &1d.

  1. Eric Winchell's code
  2. Billy Armstrong's code
  3. Steven Crisp's code
  4. Aleah Sommers code
  5. Hanna Grist's model
  6. Shahen Huda's model
  7. Kravitz' model
  8. James Stewart-Moore's model
  9. Leif Anderson's model

Modeling assignment # 2. Fault Scarp Evolution

One of the original motivations for the development of hillslope transport algorithms was in the dating of fault scarps in the American West. These researchers wished to know how long ago a fault had snapped in an earthquake, so that they could assign a level of hazards assciated with the site. Several large earthquakes had occurred in the American West that woke up the community including the Owens and Dixie Valley events:

  1. the 1872 Owens Valley quake that rattled the walls of Yosemite observed by John Muir (
  2. the 1954 Dixie Valley event in Nevada (,

In this modeling assignment, we will work a forward problem, taking the hillslope and letting it 'relax' using algorithms for the conservation of mass and for hillslope motion. You will see how these conspire to yield a diffusing form. We will treat two cases: one of a simple hillslope that begins with an abrupt offset, the other with an active fault.

Assignment 2A
Set up an initial 1D landscape z(x) in which there is a simple set of 3m height interrupting a slope with a specified gentle angle (say 4 degrees, as is typical for an alluvial fan. Let it evolve under the assumptions of conservation of colluvium, a linear relationship between colluvial discharge and local slope, and uniform material properties (i.e. the depth to bedrock is so large, so we do not have to worry about weathering it to produce mobile material). We assume a landscape diffussivity k = 0.01 m2/year.

>> Plot the evolution through time as a set of profiles stacked on one plot. (use the 'hold on' command in Matlab). Make the initial and final profiles bold.
>> Track the steepest slope at any instant and report in a plot the maximum slope versus time. See if it fits a square root of time dependence by plotting it in log-log space and inspecting the slope of the plot.

Assignment 2B
Now allow for motion on a fault that intersects the slope. Prescribe it to be a normal fault with a dip of 50 degrees, in the same direction as the dip of the surface, and allow it to slip 1m/event, with an interval of your choice. Play with this once you get it working, but we ask specifically for a uniform interval. assy with a repeat time of 1000 years (yielding a longterm slip rate of 1mm/year).

>> Again, plot profiles stacked on one plot.
>> Plot slope as a function of time, with enugh points to allow resolution of interseismic relaxation of the scarp (i.e. intervals of well less than 1000 years).
>> Make a movie of one of your simulations so that we can watch the landscape evolve.

Answers to Fault Scarp Evolution assignment # 2B

Andy's Introduction_to_Python Page features a hillslope model

Modeling assignment # 3. Settling velocities of spheres in a fluid

Here we model the simple case of single spheres, and develop settling velocity graphs and codes from which they were generated that you can use in any settling problem in the future.
Here are some key references:

  • Dietrich, W.E., 1982, Settling velocity of natural particles: Water Resources Research 18: 1615–1626.[1]
  • Ferguson, R.I. and M. Church, 2004, A SIMPLE UNIVERSAL EQUATION FOR GRAIN SETTLING VELOCITY, Journal of Sedimentary Research74: 933–937. [2]
  • Morsi, S. A. and A. J. Alexander, 1972, An investigation of particle trajectories in two-phase flow systems, Journal of Fluid Mechanics 55: 193–20. [3]


Assignment 3A
Develop a code in which you numerically integrate the equations of motion for a specified isolated sphere settling in a fluid. For now take it to be quartz in water. Report the time evolution of the vertical speed, w(t). Report three cases, with D = 10μm, 1 mm and 10 cm. Plot as well the analytical expression for this history that we developed in class, which should work well at least for the low Re (small grain size) case. Find a piece of code in which the relationship between drag coefficient, Cd, and Reynolds number, Re, is tabulated. (from Morsi and Alexander, 1972). WILL POST THIS HERE IN A COUPLE OF DAYS (02/17/2013).

Assignment 3B
Now loop through many particle diameters to develop a plot of settling speed as a function of particle diameter for the following combinations. (You will have to look up the appropriate values of viscosity and density… In each case take D to go from 1 micron to 1 m, with 5 logarithmically spaced increments in each decade of particle size. The matlab command logspace will be helpful in setting this up.) The cases are: quartz in air rain in air snow in air (or rather graupel, so that you may still treat the snow-grains as spheres) quartz in water (at 20°C) air bubble in water >> Plot your relationships in log-log space, so that we can see whether they conform to the expected D1/2 and D2 powers relevant for high and low Re developed in class. For your own ease of reference later, you may wish to plot all of these on the same plot.

In this problem you will have to determine when the particle has reached this magical (unattainable) settling speed, and at that point go on to the next particle. Unless you want it to take forever, you should allow the time step, dt, to vary from one case (particle size) to the next. Basically the idea is to take this dt to be some small fraction of the natural response time of the particle to changes in fluid velocity. Calling this latter time scale, t*, you will want to choose dt = αt*, where α is small, say 0.05, so that after 20 time steps you have gotten to the response time or characteristic time for that combination of particle/fluid.

Assignment 3C
For the ambitious, consider the following situation and answer the question posed. You are asked to create a rainfall simulator, in which raindrops are generated at a nozzle mounted high on a scaffolding over a plot in which we wish to address issues like the transport of particles by rain in a fire-scarred area in the foothills of the Front Range. You are spending a lot of effort and expense to generate a proper pdf of raindrop sizes by buying expensive nozzles. Now you are asked by your advisor, a brilliant but somewhat eccentric woman, whether the nozzle is far enough above the ground to allow the raindrops to achieve their true settling speeds. If they do not, then the kinetic energy they impart to the ground will be lower than you are expecting, with attendant departure from the theoretical calculation of sediment transport. The raindrop sizes range from 10μm to 3mm.

Question with 3C
How high should the rain experiment scaffolding and nozzle be to assure that all raindrops reach their settling speeds?

Modeling assignment # 4. Flow Velocity and Overland Flow

4A Mississippi River Velocity Structure and Sediment Transport
We are given a set of measurements of the flow velocity within a profile of the Mississippi River, at the time flowing with an approximate flow depth of 8.5m. The measurements are:

z (m above bed) U (m/s)
7.7 1.25
5.2 1.2
3.2 1.05
1.8 1.0
0.9 0.9
0.6 0.8

Recalling the law of the wall:

[math]\displaystyle{ U=\left (u*\kappa\right)ln\left (z / z_{0}\right) }[/math] (1)

wherein [math]\displaystyle{ \kappa }[/math] is the von karman constant (=0.4).

What is your best estimate of the shear velocity u*?
What do you find for the roughness length of the flow z0?

>> Plot your results in both linear-linear and log-log fashion (linear in U as the x-axis and log of z on the y-axis), showing both your data and the predicted velocity structure, using your best estimates of u* and z0. Given that:

[math]\displaystyle{ u*=\left (gHS\right)^{\frac{1}{2}} }[/math] (2)
 What is the slope of the river in this reach? Note that this is actually a really difficult measurement to make!

4B Overland Flow

Consider a hillslope of length 100m wit a uniform slope of 0.1 (i.e. 10m/ 100m). We wish to know what flow thicknesses develop on the slope in an intense rainstorm. Ultimately we are interested in knowing where on the slope the shear stress imparted by the flow is sufficient to entrain sediment.
Let the rain fall with an intensity R=40mm/hr from t=0 to t= 1hour. for now allow the infiltration rate I to be steady at a prescribed rate of 10 mm/hr (later we may allow to relax this assumption and allow it to evolve over the course of the storm). Use the Darcy-Weisbach formulation of the mean velocity U in a uniform flow (and recall Q=UH).

Calculate the evolution of the thickness of overland flow through the rainstorm, and following the end of the storm untill all of the water has drained from the slope. 

Check the thickness profile H(x) against that you calculate form the analytical solution to teh steady problem.
Report the discharge history (per unit width of the slope) as a function of time through and after the storm.