Skip to main content

NNESMO Practical 2: Evolution of river profiles

NNESMO Practical 2: Evolution of river profiles

Dr Lara Kalnins

School of Geosciences

University of Edinburgh

Contents

The changing elevation along a riverbed is referred to as the ‘longitudinal river profile’. We are interested in how this profile evolves, particularly in response to landslides creating a new cliff in the profile or to gradual landward erosion by wave action. We are also interested in how the evolution depends on factors such as the overall erosion rate, variation in catchment area, and how vulnerable steep areas are to increased erosion rates. This process can be modelled by the equation

where is distance from the head of the river, is elevation, is time, and , and are different constants that reflect different factors that influence erosion patterns.

Practical 2, Part 1: Modelling prepration

  • 2.1 – Develop a plan for how to model this system. As it is a new situation to model, it would probably be helpful to sketch the system and annotate your sketch with key dimensions, parameters you will need to know, etc.
  • 2.2 – What do you think you might use as an initial condition?
  • 2.3 – As we are working with a new equation, spend some time thinking about how the different influences mentioned in the introductory paragraph are related to the erosion equation.
  • 2.4 – Looking at the equation, why is it important that is at the head of the river rather than at the mouth? Similarly, why is it important that we take the absolute value of the gradient? Hint: consider different possible valuies of .

Practical 2, Part 2: An initial model

  • 2.5 – A reasonable test model would be a river 20 km long with an elevation drop of 1000 m. What range of values do you think would give you a reasonable resolution for this problem?
  • 2.6 – For how long do you think you are likely to need to run your model? How will this influence your choice above? Additionally, what might indicate that you are running your model for too long or too short a time for the model to be appropriate?
  • 2.7 – Now open a new script and copy in the following code, filling in the values for the different parameters and unit conversion. Sensible starting values for the constants are (rate per Myr, precise units depend on ), (unitless), and (unitless). You should also set up vectors of distances and times using and .
function NESMORiverProfile

%Put a header comment here explaining what the code does and how to use it
%%%set parameters
tmin =
tmax =
dt =
xmin =
xmax =
dx =
z_mouth =
z_head =
rate =
drainage =
power =
%%%Set up work
%unit conversions
timeconvert =
tmin_s = tmin*timeconvert;
tmax_s = tmax*timeconvert;
dt_s = dt*timeconvert;
rate_s = rate/timeconvert;
%set up arrays of distances and times, work out number of t and x steps
%needed
dist = xmin:dx:xmax;
time = tmin_s:dt_s:tmax_s;
n_end = numel(time); %n is the loop index variable for t
  • 2.8 – The simplest initial condition is a linear slope. Use the following code template to implement this, making sure it will adjust properly if you change z_min, x_max, etc. Make a plot of the initial condition.
%%%initial conditions
%Linear
%add code to set up initial condition
figure(1)
clf
subplot(2,2,1)
%add code to plot data
hold on
xlabel('Distance (m)')
ylabel('Elevation (m)')
title('Evolution of Linear Profile')
  • 2.9 – Discretise the erosion equation given above and write a subfunction to implement it. For stability, use a forward-looking discretisation across one . A central discretisation across (one in front and one behind) would in principle be more accurate (second order accurate instead of first order accurate), but the forward-looking one is much easier to work with in this equation. You will learn more about this in Session 6. Make sure your subfunction has a header comment explaining what the subfunction does and how to use it, especially what each argument does.
  • 2.10 – Add a loop over time to the main body of your code, and call your subfunction to calculate the erosion at each timestep. Using a different colour, add the results of the erosion to the plot of the initial condition. Use the following two pieces of code to make it so that the code automatically plots another profile approximately every years, where is easily set by the user. Explore the influence of changing each of the three different erosion constants.
%plotting parameters; set at beginning of script
plot_freq = %How often to plot, in years
mod_num = floor(plot_freq/dt)%Code to make the figure plot every mod_num timesteps; see plotting in main loop
%plot inside timeloop
if not((mod(n,mod_num))) %n is the index variable of the time loop
 figure(fig_num)
 subplot(2,2,1)
 %add code to plot data
end

Practical 2, Part 3: More complex initial conditions

  • 2.11 – Suppose instead of a linear initial condition, we want to model a situation where a large landslide near the mouth of the river has created a steep cliff with a (near) horizontal plane in front of it. (Landslides often have this ‘scoop’ shape, with a steep headscarp and than an essentially flat section.) Add a second model to your code where the original linear profile is modified by a steep line down to z_min, representing the headscarp, and then stays at z_min the rest of the way to the mouth. Let the user control the position of the top and base of the scarp.
  • 2.12 – Add another set of subplots showing the evolution of this second situation.
  • 2.13 – This model can generate a wider range of erosion patterns than the simple linear model. See if you can create (1) the classic exponentially flattening ‘mature’ river profile, and (2) a migrating ‘knickpoint’, where the steep section moves upstream much more rapidly than the river overall is downcutting. What sorts of parameter values create these patterns? What do you think that might say about the geological processes behind them?
  • 2.14 ADVANCED – Another common process that affects the baselevel of a river is gradual landward erosion by wave action, creating a flat, wave-cut platform where once there was a slope. Add a third model to your code that incorporates this effect, again with its own set of figures/subfigures. A reasonable rate of cliff retreat is 1 mm/yr.
  • 2.15 ADVANCED – Experiment with an initial condition of your own design.

Practical 2, Part 4: Fitting Observations

  • 2.16 – Use xlsread to load the data for the Makaha River on the island of Kauai (Seidl et al., 1994). Add a new figure showing both the river elevation and the nearby palaeosurface. Because these data are based on discrete measurements, use a symbol such as ‘x’ or ‘o’ in addition to a continuous line.
  • 2.17 – Use interp1 to create evenly spaced points along the 5 Myr old palaeosurface. We don’t use the raw palaeosurface data because it is numerically easier if the initial condition has equally spaced points, so that is constant. The data typically end a bit above sea-level, so you will need to add a point for the base of the scarp and another for the boundary condition.
  • 2.18 – Experiment with different values for , , and to try to match the present day elevation of the Makaha River.
  • 2.19 ADVANCED – Repeat your analysis for the Kauhao, Hikimoe, and Kaulaula Rivers. How similar are the values of , , and for the four rivers?

Practical 2, Part 5: Implicit timestepping – Backward Euler

In session 1, we saw that when we use the Forward Euler method, this often limits the size of , or the model becomes unstable. If , then , and we must choose a less than this for our model to be fully stable. For , the critical timestep is smaller. For models that need to run over a very long time or models that require a small to resolve the features of interest, this becomes a substantial limiting factor when using explicit timestepping methods like the Forward Euler. Instead, it is common to use implicit methods, such as Backward Euler or Crank-Nicholson.

  • 2.20 – Verify the critical timestep for your Forward Euler model. How sensitive is the critical timestep to ? To ?
  • 2.21 – Using the lecture notes from Sessions 1 and 2, discretise the power law erosion using the Backward Euler method. For the case , convert this to a set of simultaneous equation, and then express that system in matrix-vector notation. (If , then the equation is nonlinear and cannot be expressed in the form .)
  • 2.22 – Make a copy of your Forward Euler code to modify for the Backward Euler case. In the set up section, add the code to create the matrix . You will find the command diag useful, and do not forget to implement the boundary conditions.
  • 2.23 – Modify the timestepping part of your code to use the Backward Euler method. The code will no longer become unstable, but how big a timestep can you use before the answer becomes noticeably rough and inaccurate?

References

  • Seidl, M. A., Dietrich, W. E., Kirchner, J. W., 1994. Longitudinal profile development into bedrock: An analysis of Hawaiian channels. The Journal of Geology 102 (4), 457-474.
  • Stock, J. D., Montgomery, D. R., 1999. Geologic constraints on bedrock river incision using the stream power law. Journal of Geophysical Research: Solid Earth (1978-2012) 104 (B3), 4983-4993.
  • van der Beek, P., Bishop, P., 2003. Cenozoic river profile development in the Upper Lachlan catchment (SE Australia) as a test of quantitative uvial incision models. Journal of Geophysical Research: Solid Earth (1978-2012) 108 (B6).
  • Whipple, K. X., 2004. Bedrock rivers and the geomorphology of active orogens. Annu. Rev. Earth Planet. Sci. 32, 151-185.
  • Whipple, K. X., Tucker, G. E., 1999. Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs. Journal of Geophysical Research: Solid Earth (1978-2012) 104 (B8), 17661-17674.
  • Whittaker, A. C., Cowie, P. A., Attal, M., Tucker, G. E., Roberts, G. P., 2007. Bedrock channel adjustment to tectonic forcing: Implications for predicting river incision rates. Geology 35 (2),103-106.


Published with MATLAB® 7.13