Category Archives: dissertation research

Building a simple delta numerical model: Part I

This is the first in a series I’m going to write over an undetermined-number-of-parts series where we will develop all the pieces to a simple delta numerical model. I’ll outline all the pieces here while I work up an interactive version with a GUI in Matlab, and we write pieces of the code together within the posts.



The water depth (H) over a channel bed with a fixed-elevation downstream boundary can be calculated by a mass-and-momentum conservative derivation following from the full Navier-Stokes equations (e.g., Chow 1959). The so-called backwater equation, in one form, is presented as (Equation 1)
backwater_fixed_eqnwhere H is the flow depth, x is the streamwise coordinate of the flow, S is the channel bed slope, Cf  is the dimensionless coefficient of friction (corresponding to form drag in the channel bed and walls (i.e., roughness and bedforms)), and Fr is the Froude number, which scales with the ratio of inertia of the flow to the gravitational field driving the flow and can be given as (Equation 2)where qw = Qw is the width averaged water discharge, and g is the gravitational constant. Equations 1 and 2 are derived for a width-averaged, constant flow width scenario.

Equation 1 is useful to the research we do as scientists of lowland depositional systems because it predicts changes in flow depth through a delta system. The equation is solved with a boundary condition of a known flow depth (H) at the fixed-elevation downstream “receiving basin” boundary. For variations in river conditions (i.e., river water discharge), the surface elevation of receiving basins of large lowland rivers are effectively fixed, therefore the equation can be used to estimate a water surface over the lower reaches of the rivers we study for any input discharge conditions!


With a calculation of H, we can then use a conservation of water equation Qw = U A (where U is flow velocity and A = B H is channel cross sectional area for an approximately rectangular channel), to evaluate the flow velocity at all points along the channel. Calculation of U in turn allows for calculation of sediment transport, and ultimately sediment deposition — but more on these parts in future parts of this series!

Equation 1 is valid for a solution where width (B) is fixed over the downstream length of x. Where width is roughly constant, say within the lower reaches of a channel we can use this simple form of the backwater equation. However, where width is not constant over all x (i.e., dB/dx != 0), say in the river plume beyond the channel mouth, we must retain an addition term from the derivation from the Navier-Stokes equations to account for this change in width. Since we’re using computers to solve the equations, and the dB/dx term is basically only one extra gradient calculation, this will an easy addition to our model in the future. I’ll write a separate post and link it here with a solution for that version of the backwater equation as well.

So how do we solve Equation 1? Let’s see what we know: the elevation of the receiving basin is fixed at z = 0. Therefore, if we know the channel bed profile, we know at the downstream boundary for ALL possible values of Qw. We can evaluate Equation 1 from the known downstream boundary, and knowing all variables on the right hand side at this boundary, calculate the change in flow depth moving upstream over the channel bed (dH/dx)!

The solution to Equation 1 (or in its slightly more complex form for a varying width domain) can be numerically estimated pretty easily by a predictor-corrector scheme. We’ll work up a calculation for a fixed-width form of the backwater equation (Equation 1) in the next post, and then proceed to use this as one module of our simple delta model.

Building a simple delta numerical model: Part II


This material is based upon work supported by the National Science Foundation (NSF) Graduate Research Fellowship under Grant No.145068 and NSF EAR-1427177. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Dietrich Settling Velocity Matlab code

Just about two years ago, I published a code for calculating the settling velocity of a particle in water by the Dietrich, 1982 method. I later realized an error in the code that made it produce erroneous estimates. I replaced it with a corrected version of the code written in R, that you can still find on my Github page, but won’t link here. I’ve taken down the old Matlab code and unpublished the original post.

Since the first publication two years ago, I’ve dramatically improved my technical coding skills, and, just as importantly, learned a ton about how to share code with others. A new version of the code can be found linked a the end of this post, but first some philosophy about coding.

As scientists, we want to be able to reliably reproduce results and we have an attraction to openness and access to knowledge, as such, we inevitably share code with one another. I’ve shared a good bit of my code thus far in my career, and I’ve made a few notes about the matter below:

  1. It is better to organize your code in a way that may not be the most succinct (when speed is not an issue) while improving the readability of the code. In essence, making your code readable to those unfamiliar with it, so they can follow along with your logic.
  2. Comment everything. It is so easy to get into the groove of writing something and blow through it until it’s finished and then never return to clean up and document the code. When others (or even you) go back to look at it, simple comments explaining what a line does go miles further than a detailed explanation in the documentation.
  3. Write your code to withstand more that you are currently using it for. What I mean by this is to not write code that works only when you are manipulating this particular dataset or working this particular problem, but code that will work over a broader range of data or can be applied to another problem. Along the same line, use explicit values as sparingly as possible, e.g., if X = [1:12] is a 12 element vector and you want the middle value, do not use X(6), but instead use X(length(X)/2); or slightly better yet, X(ceil(length(X)/2)) to at least not throw an error and give you the first of the two middle rows if X is an odd number of elements in length.


Anyway, the code is uploaded to my Github and you can find it there. The documentation output on a help get_DSV call in Matlab is produced below.



Dietrich, W. E. Settling velocity of natural particles. Water Resources Research 18, 1615–1626 (1982).

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No.1450681. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation.

Yellow River shoreline traces

For some ongoing research I needed a set of shoreline traces of the Yellow River delta that I could measure some properties of to compare modeling results to. My model predicts delta lobe growth rates, for example the growth of the Qingshuiguo lobe on the eastern side of the delta from 1976 to 1996.

early growth stage of Qingshuiguo lobe of the Yellow River delta, figure from van Gelder et al., 1994.

early growth stage of Qingshuiguo lobe of the Yellow River delta, figure from van Gelder et al., 1994.

It’s a little early for me to release the data I’ve collected from the shoreline traces (n=296), or the code I used to generate the traces — but below is a figure mapping all the shoreline traces I have from 1855 to present. My model can also predict rates of delta system growth (i.e., growth over many lobe cycles) and I have generated data from these traces that I can compare to as well. I like the below figure, because I think it nicely demonstrates the rapid reworking of delta lobes (yes even on “river dominated” deltas!); evidenced most clearly by the retreat of the Diakou lobe to the north side of the delta.



This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No.1450681. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation.

Yellow River (Huanghe) 2016 survey plans

I’ll be departing soon for my second field campaign on the Yellow River, China. We traveled there around the same time last year to generate a dataset we could begin to explore and develop research ideas; now, a year later, we will return with a more focused plan to gather the data we need to address our research questions. We will bring with us a suite of sophisticated data collection equipment. Below is a terrible sketch I made to demonstrate most of the data we will collect at each of our stations located on the river.


Pretend that the grey thing floating on the water might look like a boat, viewed from the back. We aim to characterize the mechanisms for sediment transport within the Yellow River at various stages of a flood discharge curve. We use a point-integrating sediment sampler to collect numerous suspended sediment and water samples from various depths in the water column (yellow stars). In order to compare the compositions of suspended material to its source (i.e., the bed), we use a Shipek grab sampling device to collect sediment from the bed for our analysis. Deployed off the other side of the boat, we will can characterize the velocity structure of the water column through the use of an aDcp (acoustic velocity profiler) and a mechanical propeller driven velocimeter for near-bed measurements where the aDcp may lack resolution.

This setup represents only a portion of our survey plans, of course we’ll be there for six weeks that will generate a diverse and (hopefully!) comprehensive dataset for our future work.

GSA abstract acceptance!

A little while ago, I heard back from the annual conference planners for Geological Society of America that my abstract submission was accepted for presentation! I’ll be presenting the following poster at GSA 2015 in Baltimore in early November.

Paper No. 311-3
Presentation Time: 9:00 AM-6:30 PM


MOODIE, Andrew J.1, MA, Hongbo1, NITTROUER, Jeffrey A.2, CARLSON, Brandee N.3 and KINEKE, Gail C.4, (1)Earth Science, Rice University, 6100 Main Street, MS-126, Houston, TX 77005, (2)Dept of Earth Science, Rice University, 6100 Main Street, MS-126, Houston, TX 77005, (3)Department of Earth Science, Rice University, 6100 Main Street MS-126, Houston, TX 77005, (4)Earth and Environmental Science, Boston College, 140 Commonwealth Ave, Chestnut Hill, MA 02467

The Huanghe (Yellow River), China, is an end-member fluvial system because it is characterized by very high sediment loads, variable water discharge over annual and decadal timescales, and substantial anthropogenic influences. Over the past several centuries, the channel bed of the Huanghe has aggraded significantly, resulting in several catastrophic avulsions. In order to mitigate channel-bed aggradation, in 2002, Chinese engineers implemented an annual Water and Sediment Discharge Regulation (WSDR) period when, over a two-week period, a controlled flood wave designed to erode the channel bed is triggered by releasing water from a network of upstream dams and reservoirs. Here, we present the results of an extensive time-series field survey that measured fluid and sediment properties to assess the effect of a WSDR flood wave on a portion of the river located ~80 km upstream of the ocean outlet. This controlled experiment represents a classic external perturbation to a fluvial-deltaic system; the morphodynamic responses, however, were anything but predictable. For example, multibeam bathymetric data illuminate the uniqueness of channel-bed properties, including a broad, flat, and shallow sand-bed devoid of a thalweg or bedforms in the straight-reach segments of the Huanghe. Alternatively, in bend segments, the channel-bed deepens significantly, and longitudinal dunes emerge as the dominant bed features. Interestingly, as the flood wave progresses over time, the channel bed in bend segments rapidly erodes upstream and downstream, thus propagating dunes in both directions. In order to further ascertain linkages between the fluid flow, sediment transport, and morphological evolution of the bed, we couple water velocity, suspended sediment concentration, and bedload transport datasets to standard physical models. The results are then used to constrain a numerical model that assesses stability of the channel bed over the duration of the flood wave. In effect, this study shows the inherent unpredictability that emerges between contemporaneously evolving “external” (allogenic) and “internal” (autogenic) processes. Finally, the results of this field and numerical investigation are used to constrain a prediction for the long-term stratigraphic development of the Huanghe fluvial-deltaic system.

One year of graduate school travel

The 24th of August 2015 marks my one year anniversary of beginning graduate school. It has been an awesome year, filled with new science, great friends, and lots of travel. I kept a record of where I slept every night for the past year to see just how far I made it. I visited three continents (plus Central America), six states, and 32 different cities. I spent just over 30% of the year away from Houston (112 nights). I got to see some really amazing places on the planet. It puts into perspective for me how just how great of a gig being a grad student is — get paid to do cool science.

Below, I made some maps to show where I went. The color of a dot represents when I traveled there during the year (Aug 24 2014 to Aug 24 2015) according to traditional rainbow color order. The maps were made in GMT using this script.


Blow ups:
global map
US map
Gulf Coast map

Cheers to what next year brings!

Back from China!

I’ve finally returned from China after a 6-week field campaign to get a first round look at my field location. We worked hard every single day after our gear was released by China Customs officials on the Yellow River Delta doing Vibracoring for Brandee Carlson’s research and doing Multibeam and river survey work for my work. Below you can see me running our Multibeam setup aboard our research vessel. More pictures and descriptions to come!


River hysteresis and plotting hysteresis data in R

Hysteresis is the concept that a system (be it mechanical, numerical, or whatever else) is dependent of the history of the system, and not only the present conditions. This is the case for rivers. For example, consider the following theoretical flood curve and accompanied sediment discharge curve (Fig. 1a).

Figure 1. Theoretical plots to demonstrate the hysteresis of a river.

Figure 1. Theoretical plots to demonstrate the hysteresis of a river.

With the onset of the flood, the increased sediment transport capacity of the system entrains more sediment and the sediment discharge curve (red, Fig. 1a) rises. However, the system may soon run out of sediment to transport (really just a reduction in easily transportable sediment), and the sediment discharge curve decreases although the water discharge curve remains high in flood.

In Fig. 1b, the sediment discharge and water discharge are plotted through time, a typical way of observing the hysteresis of a system. Note that for the rising limb and falling limb of the river flood, the same water discharge produces two different sediment transport values.

Now, let’s imagine that we want to investigate how important the history of the system is to the present state of our study river. You can grab the data I’ll use, here. This is data from one year of flow on the Huanghe (Yellow River) in China, and it has been smoothed with a moving average function to make the hysteresis function more visible.

Making the plot

It is easy enough to plot a line with R (the lines function) but with a hysteresis plot, it is important to be able to determine which direction time is moving forward along the curve. For this reason we want to use arrows. So we plot the line data first, with:


and then using a constructed vector of every 22nd number, we plot an arrow over top of the lines using:

s <- seq(from=1, to=length(df[,"Qw"])-1, by=22)
arrows(df[,"Qw"][s], df[,"Qs"][s], df[,"Qw"][s+1], df[,"Qs"][s+1], 

Finally, with a few more modifications to the plot (see the full code here), we can produce Fig. 2 below. This plot is comparable to the theoretical one above.

Figure 2: Hysteresis plot from actual data.

Figure 2. Hysteresis plot from actual data.

Using the green lines and points, I have highlighted the observation that for the rising limb and falling limb of a flood, there can be substantially different sediment discharges for the same water discharge — this observation is not so easily made from the plot on the left, but it is immediately clear in the hysteresis plot on the right.