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.

5 thoughts on “River hysteresis and plotting hysteresis data in R”

  1. Dear Andrew, thank you very much for the detailed instructions!

    Can you please explaine a bit more about smoothing with a moving average function? How did it changed the original data?

  2. Hi Andrew, I am a graduate student trying to use R to calculate the Hysteresis Index (HI) proposed by Lawler et al. 2016 https://www.sciencedirect.com/science/article/pii/S0048969715310093
    Thank you for showing how to plot the hysteresis loops, is there a way to calculate the Hysteresis Index in R by subtracting the falling limb of the curve from the rising limb for each 5% percentile of the data and summing them? I am trying to compute the HI values for my normalized water discharge and turbidity data, thanks!

    1. Hi Junjie, I have not written a function to do what you describe. Though it sounds like it should be possible, and reasonably simple. The approach I would take would be to 1) split the data to rising and falling data, 2) sort the values by magnitude, 3) compute the CDF of the data, loop through each 5% of the CDF and compare between rising and falling. Good luck!

    2. Anatolii Tsyplenkov

      Dear Junjie Chen,

      I’ve coded HImid index from [Lawler et al., 2006] previously. As I understand you can just change 0.5 to 0.05 in the second raw, i.e. Qmiddle discharge to Q0.05

      himid <- function(x){

      mid <- function(x) 0.5*(max(x, na.rm = T) – min(x, na.rm = T)) + min(x, na.rm = T)

      target <- mid(x$q)
      idx 0)

      f <- function(i, target) approx(c(x$q[i], x$q[i+1]), c(x$ssc[i], x$ssc[i+1]), xout=target)$y
      yp <- sapply(idx, f, target = target)

      if (!length(yp) == 2) {
      h <- NA
      } else {

      if (yp[1] < yp[2]) {
      h <- -1/(yp[1]/yp[2]) + 1
      } else {

      h <- (yp[1]/yp[2]) – 1


Leave a Reply

Your email address will not be published. Required fields are marked *