4  Stage-Discharge Rating Curves

The power function for relating discharge and water level is:

\[ Q = K(H - H_0)^z \]

where:

\(Q\) = discharge in cfs; \(H\) = stage in feet; \(H_0\) = stage at zero flow (in feet); \(K\) and \(z\) are parameters determined by fitting stage-discharge measurements.

Sometimes we know \(H_0\) based on observed data and can plug the known value in. It is also possible to estimate it using regression. Here I am going to let the non linear least squares solve for the best value because we are likely going to fit multiple models anyways. You should check to make sure that the the solved parameter makes sense (not negative or highly positive). \(K\) is equal to the discharge at which \((H-H_0) = 0\) and \(Z\) describes the slope. Using the nls() function requires the analyst provide some starting parameters, the choice of starting parameters is important because the model can converge on the wrong solution or fail to find a solution if the parameters are far away from the optimal answer. Using the nls_multstart() function from the nls.mutlstart package, we can provide a range of possible starting values and the function will iterate combinations of starting values and select the best model using information theory metrics.

4.1 Data

This example uses field measurements of stream stage and instantaneous discharge measured at the COMO NEON field site (NEON (National Ecological Observatory Network), n.d.). Raw data is available in the neon.rds file within the tutorial download.

library(tidyverse)
library(twriTemplates)

## neon_stage_discharge is an example dataset in the twriTemplates package
stage_discharge <- neon_stage_discharge
glimpse(stage_discharge)
Rows: 134
Columns: 36
$ uid                     <chr> "1068e9ad-37bd-42ac-96d2-b999545adefe", "06d63…
$ recorduid               <chr> "3b2d4830-50cd-40c8-a772-790d2386c33d", "5ed72…
$ domainID                <chr> "D13", "D13", "D13", "D13", "D13", "D13", "D13…
$ siteID                  <chr> "COMO", "COMO", "COMO", "COMO", "COMO", "COMO"…
$ namedLocation           <chr> "COMO.AOS.discharge", "COMO.AOS.discharge", "C…
$ collectedBy             <chr> "LCLARK", "HSCHARTEL", "DMONAHAN", "HSCHARTEL"…
$ startDate               <dttm> 2015-08-11 15:30:00, 2015-08-24 16:16:00, 201…
$ collectDate             <dttm> 2015-08-11 15:30:00, 2015-08-24 16:16:00, 201…
$ streamStage             <dbl> 0.30, 0.34, 0.33, 0.31, 0.32, 0.35, 0.37, 0.74…
$ streamStageUnits        <chr> "meter", "meter", "meter", "meter", "meter", "…
$ handheldDeviceID        <chr> "122581001261", "122581001261", "122581001261"…
$ velocitySensorID        <chr> "132660300436", "132660300436", "132660300436"…
$ filterParamTime         <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 20, 10, 10…
$ stationEntryTest        <chr> "Non-fixed", "Non-fixed", "Non-fixed", "Non-fi…
$ flowCalculation         <chr> "Mid-section", "Mid-section", "Mid-section", "…
$ waterEdge               <chr> "Right edge water", "Right edge water", "Right…
$ totalDischarge          <dbl> 10.81, 5.71, 6.58, 3.49, 3.15, 9.39, 12.75, 77…
$ totalDischargeUnits     <chr> "litersPerSecond", "litersPerSecond", "litersP…
$ samplingProtocolVersion <chr> "NEON.DOC.001085vB", "NEON.DOC.001085vB", "NEO…
$ averageVelocityUnits    <chr> "meterPerSecond", "meterPerSecond", "meterPerS…
$ waterDepthUnits         <chr> "meter", "meter", "meter", "meter", "meter", "…
$ tapeDistanceUnits       <chr> "meter", "meter", "meter", "meter", "meter", "…
$ flowCalcQF              <chr> "0", "0", "0", "0", "0", "0", "0", "0", "0", "…
$ dischargeUnitsQF        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ streamStageUnitsQF      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ averageVelocityUnitsQF  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ waterDepthUnitsQF       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ tapeDistanceUnitsQF     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ lowVelocityFinalQF      <dbl> 5, 20, 0, 10, 0, 17, 0, 45, 0, 11, 25, 11, 6, …
$ finalDischarge          <dbl> 10.71, 5.62, 6.52, 3.44, 3.07, 9.29, 12.75, 78…
$ totalDischargeCalcQF    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ profileName             <chr> "COMO", "COMO", "COMO", "COMO", "COMO", "COMO"…
$ stageImpractical        <chr> "OK", "OK", "OK", "OK", "OK", "OK", "OK", "OK"…
$ dataQF                  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ publicationDate         <chr> "20211222T003739Z", "20211222T003739Z", "20211…
$ release                 <chr> "RELEASE-2022", "RELEASE-2022", "RELEASE-2022"…

4.2 Plot Data

First step is to plot the data.

# plot observations over time
ggplot(stage_discharge) +
  geom_point(aes(collectDate, finalDischarge)) +
  labs(x = "Date", y = "Discharge [L/s]") +
  theme_TWRI_print()
Warning: Removed 9 rows containing missing values (`geom_point()`).
# plot the stage and discharge relationship
ggplot(stage_discharge) +
  geom_point(aes(streamStage, finalDischarge)) +
  labs(x = "Stage [m]", y = "Discharge [L/s]") +
  scale_x_log10() +
  scale_y_log10() +
  theme_TWRI_print()
Warning: Removed 9 rows containing missing values (`geom_point()`).

Looks like there might be a shift in the rating curve at some point. Use dplyr to create a new year variable and the gghighlight package to explore it.

# install.packages("gghighlight")
library(gghighlight)
stage_discharge |> 
  filter(!is.na(finalDischarge)) |> 
  mutate(year = format(collectDate, format = "%Y")) |> 
  ggplot() +
  geom_point(aes(streamStage, finalDischarge, color = year)) +
  facet_wrap(~year) +
  gghighlight(max_highlight = 8,
              use_direct_label = FALSE,
              calculate_per_facet = FALSE)   +
  labs(x = "Stage [m]", y = "Discharge [L/s]") +
  theme_TWRI_print()

4.3 Fit Rating Curve

It appears the rating curve shifted after 2016 and possibly after 2018. In practice we would fit different rating curves based on this information. For this example the rating curve will only be fit to 2019-2021 data. We will use this data to estimate the \(K\), \(H_0\) and \(Z\) parameters in the power function described at the top of the chapter using nonlinear least squares.

# install.packages("nls.multstart")
# install.packages("nlstools")
library(nls.multstart)
library(nlstools)

# clean the data a little bit and filter
stage_discharge <- stage_discharge |> 
  filter(!is.na(finalDischarge)) |>
  mutate(year = format(startDate, "%Y")) |> 
  filter(year >= 2018) |> 
  # convert units to feet and cfs
  mutate(streamStage = streamStage * 3.28084,
         finalDischarge = finalDischarge * 0.0353147)

# Set the equation
f_Q <- formula(finalDischarge ~ K*(streamStage - H_0)^Z)

# Some initial starting values
start_lower <- c(K = -10, Z = -10, H_0 = 0.02)
start_upper <- c(K = 10, Z = 10, H_0 = 1)

# nonlinear least squares
m1 <- nls_multstart(f_Q,
          data = stage_discharge,
          iter = 1000,
          start_lower = start_lower,
          start_upper = start_upper,
          supp_errors = 'Y',
          lower = c(K = -10, Z = -10, H_0 = 0),
          control = minpack.lm::nls.lm.control(maxiter = 1000L))

summary(m1)

Formula: finalDischarge ~ K * (streamStage - H_0)^Z

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K     5.0398     1.3422   3.755 0.000335 ***
H_0   0.7874     0.1315   5.990 6.31e-08 ***
Z     1.9098     0.2568   7.435 1.23e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.889 on 77 degrees of freedom

Number of iterations to convergence: 42 
Achieved convergence tolerance: 1.49e-08

NLS estimated parameters are: \(K =\) 5.0397806, \(H_0 =\) 0.7874016, and \(Z =\) 0.7874016. Before using these parameter, evaluated the goodness of fit using the model residuals (add citation or note with more resources here).

Show the code
# for easy multipanel plots, use the patchwork package
#install.packages("patchwork")
library(patchwork)

std_resids <- as_tibble(nlsResiduals(m1)$resi2)

stage_discharge <- stage_discharge |> 
    mutate(fits = std_resids$`Fitted values`,
         residuals = std_resids$`Standardized residuals`)


p1 <- ggplot(stage_discharge) +
  geom_density(aes(residuals)) +
  labs(x = "Standardized Residuals",
       y = "Count",
       subtitle = "Distribution of standardized residuals") +
  theme_TWRI_print()


p2 <- ggplot(stage_discharge) +
  geom_point(aes(streamStage, residuals), color = "steelblue", alpha = 0.4) +
  labs(x = "Stream Height [ft]",
       y = "Standardized Residuals",
       subtitle = "Residuals against stream height") +
  theme_TWRI_print()

p3 <- ggplot(stage_discharge) +
  geom_point(aes(fits, finalDischarge), color = "steelblue", alpha = 0.4) +
  labs(x = "Model Fits",
       y = "Measured Discharge [cfs]",
       subtitle = "Measured against fitted") +
  theme_TWRI_print()

p4 <- ggplot(stage_discharge) +
  stat_qq(aes(sample = residuals), color = "steelblue", alpha = 0.4) +
  stat_qq_line(aes(sample = residuals)) +
  labs(x = "Theoretical",
       y = "Standardized Residuals",
       subtitle = "Sample Quantile against theoretical quantile") +
  theme_TWRI_print()


# patchwork allows us to assemble plots using + and /

(p1 + p2) / (p3 + p4)

The plots indicate increasing residual variance as stream height increases and heavy tails in Q-Q plot. Two options from here are to (1) fit a rating curve per year or (2) fit a piece-wise log-linear model that assumes a different relationship along different parts of the rating curve. Because of the large residual variance at the top of curve, I think fitting a curve per year will do the job. You may also choose to fit models by season (small streams in particular may see large changes in rating curves do to changes in stream bank vegetation).

4.4 Fit Multiple Curves

Note

This section is slightly more advanced and requires some understanding of list structures in R and nested dataframes.

The purrr package facilitates running a function on a list of nested data. The idea here is to subset the dataframe by year, create a list with each item in the list being a subset of the dataframe, then running nls_multstart() on each item in that list. Sounds like a loop huh? We achieve this using for() or lapply() functions in base R. The nice thing about doing it with purrr is that we can keep everything together in a single dataframe. The first step is to create a nested dataframe:

nested_data <- neon_stage_discharge |> 
  filter(!is.na(finalDischarge)) |>
  mutate(year = format(startDate, "%Y")) |> 
  filter(year >= 2018) |> 
  mutate(streamStage = streamStage * 3.28084,
         finalDischarge = finalDischarge * 0.0353147) |> 
  ## group data by year
  group_by(year) |> 
  ## nest the data by year
  nest()

nested_data
# A tibble: 4 × 2
# Groups:   year [4]
  year  data              
  <chr> <list>            
1 2018  <tibble [20 × 36]>
2 2019  <tibble [24 × 36]>
3 2020  <tibble [15 × 36]>
4 2021  <tibble [21 × 36]>

Now we use the map() function in purrr to iterate the nls.multstart() function on each nested dataframe:

nested_data <- nested_data |> 
  mutate(model_output = map(.x = data,
                            ~nls_multstart(formula = f_Q,
                                           data = .x,
                                           iter = 1000,
                                           start_lower = start_lower,
                                           start_upper = start_upper,
                                           supp_errors = 'Y',
                                           lower = c(K = -10, Z = -10, H_0 = 0),
                                           control = minpack.lm::nls.lm.control(maxiter = 1000L))))
nested_data
# A tibble: 4 × 3
# Groups:   year [4]
  year  data               model_output
  <chr> <list>             <list>      
1 2018  <tibble [20 × 36]> <nls>       
2 2019  <tibble [24 × 36]> <nls>       
3 2020  <tibble [15 × 36]> <nls>       
4 2021  <tibble [21 × 36]> <nls>

We created a new column called model_output that is a list of the output from the nls_multstart() function that was run on each item in the data column. Grab the residuals and fits from each model:

Show the code
nested_data <- nested_data |> 
  mutate(residuals = map(model_output,
                         ~as_tibble(nlsResiduals(.x)$resi2))) |> 
  unnest(c(data,residuals))

p1 <- ggplot(nested_data) +
  geom_point(aes(`Fitted values`, `Standardized residuals`)) +
  facet_wrap(~year) +
  labs(x = "Fitted",
       y = "Standardized Residuals",
       subtitle = "Sample Quantile against theoretical quantile") +
  theme_TWRI_print()

p2 <- ggplot(nested_data) +
  geom_density(aes(`Standardized residuals`)) +
  facet_wrap(~year) +
  labs(x = "Standardized Residuals",
       y = "Count",
       subtitle = "Distribution of standardized residuals") +
  theme_TWRI_print()

p3 <- ggplot(nested_data) +
  geom_point(aes(`Fitted values`, finalDischarge), color = "steelblue", alpha = 0.4) +
  facet_wrap(~year) +
  labs(x = "Model Fits",
       y = "Measured Discharge [cfs]",
       subtitle = "Measured against fitted") +
  theme_TWRI_print()

p4 <- ggplot(nested_data) +
  stat_qq(aes(sample = `Standardized residuals`), color = "steelblue", alpha = 0.4) +
  stat_qq_line(aes(sample = `Standardized residuals`)) +
  facet_wrap(~year) +
  labs(x = "Theoretical",
       y = "Standardized Residuals",
       subtitle = "Sample Quantile against theoretical quantile") +
  theme_TWRI_print()

p1 / p2 / p3 / p4

The results are a mixed bag. For the most part, the residuals are tighter to mean zero and the tails are not as heavy as the first example. We will assume the data is good enough to continue the example. In practice, I might explore the use of seasonal rating curves or piece-wise functions.

Next lets make a nice plot showing the rating curve with the observed data.

fits <- nested_data |> 
  nest(data = -c(year, model_output)) |> 
  # create a new dataframe by group
  # this includes the full range of the predictor variable
  # so we can draw a nice smooth line using predictions
  mutate(newdata = map(data,
                       ~{
                         tibble(streamStage = seq(min(nested_data$streamStage),
                                                  max(nested_data$streamStage),
                                                  length.out = 100))
                       })) |> 
  mutate(fits = map2(newdata, model_output,
                     ~{predict(.y, .x)})) |> 
  unnest(c(newdata, fits))

ggplot() +
  geom_point(data = nested_data,
             aes(streamStage, finalDischarge, color = year)) +
  geom_line(data = fits,
            aes(streamStage, fits, color = year)) +
  facet_wrap(~year) +
  gghighlight(max_highlight = 8,
              use_direct_label = FALSE,
              calculate_per_facet = FALSE)   +
  labs(x = "Stage [ft]", y = "Discharge [cfs]") +
  theme_TWRI_print()