Flood Frequency Analysis

Marc Los Huertos

EA030 · Compiled 2026-04-24

EA030 — RTricks — Activity 11

Flood Frequency Analysis: Return Intervals and the Gumbel Distribution

Prof. Marc Los Huertos  •  Pomona College  •  Environmental Analysis

← Back to RTricks

Concepts Packages Settings Download Empirical Method Gumbel Fit Log-Pearson III Split Periods Helper Functions Summary

🎯 Learning Outcomes

  1. Explain what a flood return period means probabilistically and why it is not a guarantee of recurrence.
  2. Download USGS annual peak flow records using readNWISpeak() from the dataRetrieval package.
  3. Apply the Weibull plotting-position formula to rank annual peaks and estimate empirical return periods.
  4. Fit a Gumbel (Extreme Value Type I) distribution to annual maxima using the method of moments and plot results on a Gumbel probability axis with 95% confidence intervals.
  5. Apply the Log-Pearson Type III procedure to estimate flood quantiles at standard return periods.
  6. Split a long record into two sub-periods, compare the fitted frequency curves, and interpret what a shift reveals about non-stationarity.

What this activity covers. On April 5, 2025, dozens of USGS gauges across the central United States registered their highest discharge on record. Flood frequency analysis (FFA) places events like that in a probabilistic context: given the historical record at a gauge, how rare is this peak? This activity builds the complete FFA pipeline — from raw annual peak data through empirical plotting positions, Gumbel distribution fitting, and a two-period stationarity comparison — using the dataRetrieval package and base-R graphics.

1. Background Concepts

1.1 The hydrograph and flood peaks

A hydrograph shows discharge over time during and after a storm. Its key features are:

Flood frequency analysis focuses exclusively on the annual peak — the single highest discharge recorded in each water year (October 1 – September 30).

1.2 Return period — what it actually means

A T-year flood is a discharge that has a 1/T probability of being equalled or exceeded in any single year. The 100-year flood is not a flood that happens once a century; it is a flood that has a 1% annual probability of occurrence.

⚙️ Annual exceedance probability
P(Q ≥ Q_T in any one year) = 1 / T

P(Q ≥ Q_T at least once in n years) = 1 − (1 − 1/T)^n

For a 100-year flood over a 30-year mortgage: 1 − (1 − 0.01)^30 = 26%.
That is the probability a 100-year flood affects a property during a typical mortgage term.

1.3 Stationarity and its discontents

Classical FFA assumes floods are stationary — drawn from a distribution whose parameters do not change over time. Milly et al. (2008) famously declared “stationarity is dead” given evidence that land use change, reservoir operations, and climate change all alter flood statistics. Comparing frequency curves fitted to early vs. recent sub-periods of a record is a first diagnostic for non-stationarity.

📌 Annual peak vs. daily mean discharge. Activity 10 used daily mean discharge from readNWISdv(). This activity uses instantaneous annual peak discharge from readNWISpeak(). Peak instantaneous discharge can be 2–5× larger than the daily mean for the same event because the daily mean averages the entire 24-hour period including the rising and falling limbs. Always use annual peaks — not daily means — for flood frequency analysis.

2. Packages

# install.packages(c("dataRetrieval", "tidyverse", "knitr", "kableExtra", "e1071"))

library(dataRetrieval)  # USGS NWIS data access
library(tidyverse)      # dplyr, ggplot2
library(knitr)          # tables
library(kableExtra)     # styled HTML tables
library(e1071)          # skewness() for Log-Pearson III

3. Settings

All user-defined parameters are collected here. Change only these values to analyse a different station.

# ── Station settings ──────────────────────────────────────────────────────────
mysite    <- "07061000"    # Current River at Van Buren, MO
startDate <- ""            # blank = use full period of record
endDate   <- ""

# ── Reference discharge for a recent extreme event ────────────────────────────
# April 5, 2025 peak at this station (cfs)
Qmax_ref  <- 180000

# ── Date to split the record into two sub-periods ────────────────────────────
splitdate <- "2000-01-01"

💡 Choosing a split date. The split date divides the record into “historical” and “recent” periods. A good choice is a watershed-change event: a dam completion, a major land-use shift, or a climatological breakpoint identified in the literature. For the Current River, Missouri, 2000 roughly separates a period of relative stationarity from the more volatile recent decades. For your own station, check whether a change-point analysis (e.g., the changepoint package) identifies a structural break.

4. Download Annual Peak Data

readNWISpeak() queries the USGS annual peak flow file — a curated record of the single highest instantaneous discharge in each water year, maintained by USGS hydrographers who quality-control and approve every value.

annualpeak <- readNWISpeak(mysite, startDate, endDate)

# Remove records with missing peak discharge or date
annualpeak <- annualpeak %>%
  filter(!is.na(peak_va), !is.na(peak_dt))

cat("Station       :", mysite, "\n")
## Station       : 07061000
cat("Period        :", as.character(min(annualpeak$peak_dt)),
    "to", as.character(max(annualpeak$peak_dt)), "\n")
## Period        : 1904-03-27 to 2024-07-18
cat("Annual peaks  :", nrow(annualpeak), "\n")
## Annual peaks  : 95
cat("Discharge range (cfs):", round(min(annualpeak$peak_va)),
    "to", format(max(annualpeak$peak_va), big.mark = ","), "\n")
## Discharge range (cfs): 22400 to 382,000
ggplot(annualpeak, aes(x = peak_dt, y = peak_va)) +
  geom_col(fill = "#2e7d94", alpha = 0.75, width = 200) +
  geom_hline(yintercept = Qmax_ref, linetype = "dashed",
             color = "#c0392b", linewidth = 0.9) +
  annotate("text", x = min(annualpeak$peak_dt),
           y = Qmax_ref * 1.04,
           label = paste0("Reference Q = ",
                          format(Qmax_ref, big.mark = ",")),
           color = "#c0392b", hjust = 0, size = 3.4, fontface = "italic") +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title    = paste("Annual Peak Discharge \u2014 USGS", mysite),
    subtitle = "Red dashed line = reference extreme event",
    x = NULL, y = "Peak discharge (cfs)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold",
                                 color = "#1e3a2f", size = 13),
    plot.subtitle = element_text(color = "#7a7a72", size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor   = element_blank()
  )

5. Empirical Return Periods — Weibull Plotting Positions

The simplest estimate of a return period requires no distributional assumption. Rank the annual peaks from largest to smallest, then apply the Weibull formula:

⚙️ Weibull plotting position
T = (n + 1) / m

where n = total years of record, m = rank (1 = largest value).
This is an unbiased estimator of the exceedance probability 1/T under no distributional assumption.

Q <- annualpeak$peak_va
n <- length(Q)

# Rank from largest (rank 1) to smallest
r <- n + 1 - rank(Q)
T_emp <- (n + 1) / r

emp_df <- tibble(
  year       = format(annualpeak$peak_dt, "%Y"),
  discharge  = Q,
  rank       = r,
  return_period = round(T_emp, 2),
  exceed_prob   = round(1 / T_emp, 4)
) %>%
  arrange(rank) %>%
  slice_head(n = 15)

emp_df %>%
  kbl(
    caption   = paste("Top 15 annual peaks with Weibull return periods,",
                      "USGS", mysite),
    col.names = c("Year", "Discharge (cfs)", "Rank",
                  "Return period (yr)", "Exceedance probability")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)
Top 15 annual peaks with Weibull return periods, USGS 07061000
Year Discharge (cfs) Rank Return period (yr) Exceedance probability
1916 382000 1 96.00 0.0104
1915 373000 2 48.00 0.0208
1927 369000 3 32.00 0.0312
1945 324000 4 24.00 0.0417
1982 312000 5 19.20 0.0521
1943 281000 6 16.00 0.0625
1935 266000 7 13.71 0.0729
1938 260000 8 12.00 0.0833
1907 259000 9 10.67 0.0938
1927 244000 10 9.60 0.1042
1949 236000 11 8.73 0.1146
1921 224000 12 8.00 0.1250
1904 220000 14 6.86 0.1458
1918 220000 14 6.86 0.1458
1933 220000 14 6.86 0.1458 
# Gumbel reduced variate for the x-axis
y_emp   <- -log(-log(1 - 1 / T_emp))

# Reference Ttick sequence for axis labels
Ttick <- c(1.01, 1.1, 1.5, 2, 5, 10, 20, 50, 100)
ytick <- -log(-log(1 - 1 / Ttick))

plot(
  y_emp, Q,
  xaxt = "n",
  xlab = "Return period T (years)",
  ylab = "Annual peak discharge (cfs)",
  main = paste("Empirical flood frequency \u2014 USGS", mysite),
  pch  = 21, bg = "#8aad96", col = "#1e3a2f", cex = 1.1,
  ylim = c(0, max(Q) * 1.1),
  xlim = c(min(ytick), max(ytick))
)
axis(1, at = ytick, labels = as.character(Ttick))
abline(h = Qmax_ref, lty = 2, col = "#c0392b", lwd = 1.8)
abline(v = ytick, lty = 3, col = "lightgray")
abline(h = seq(0, max(Q) * 1.1, length.out = 6), lty = 3, col = "lightgray")
legend("topleft",
       legend = c("Annual peak", paste("Reference Q =",
                                        format(Qmax_ref, big.mark = ","))),
       pch    = c(21, NA),
       lty    = c(NA, 2),
       col    = c("#1e3a2f", "#c0392b"),
       pt.bg  = c("#8aad96", NA),
       bty    = "n", cex = 0.85)

📌 Reading the Gumbel probability axis. The x-axis uses the Gumbel reduced variate y = −ln(−ln(1 − 1/T)), which linearises data that follows the Gumbel distribution. If the points fall on a straight line, the Gumbel distribution is a good fit. Systematic curvature indicates a heavier or lighter tail than Gumbel predicts — try the Log-Pearson III in that case.

6. Gumbel Distribution Fit — Method of Moments

The Gumbel (Extreme Value Type I) distribution is the most widely applied parametric model for annual flood maxima. Its CDF is:

⚙️ Gumbel CDF and frequency factor
F(x) = exp[−exp(−(x − μ)/β)]

Method-of-moments estimates:
β = s · √6 / π     μ = x̄ − 0.5772 · β

Frequency factor for return period T:
K_T = −(√6/π) · [0.5772 + ln(ln(T/(T−1)))]

Flood quantile: Q_T = x̄ + K_T · s

# ── Fit parameters ────────────────────────────────────────────────────────────
xbar <- mean(Q);  s <- sd(Q);  nQ <- length(Q)

# Return period grid for fitted line
Ttick_fine <- c(1.001, 1.01, 1.1, 1.5, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                15, 20, 25, 30, 35, 40, 50, 60, 70, 80, 90, 100)
ytick_fine <- -log(-log(1 - 1 / Ttick_fine))

# Frequency factor and fitted quantiles
K_T    <- -(sqrt(6) / pi) * (0.5772 + log(log(Ttick_fine / (Ttick_fine - 1))))
Q_T    <- xbar + K_T * s

# 95% CI using Kite (1949) standard error formula
se_T   <- (s * sqrt(1 + 1.14 * K_T + 1.1 * K_T^2)) / sqrt(nQ)
LB     <- Q_T - qt(0.975, nQ - 1) * se_T
UB     <- Q_T + qt(0.975, nQ - 1) * se_T

# Axis label positions
Ttick_lab <- c(1.01, 1.1, 1.5, 2, 5, 10, 20, 50, 100)
ytick_lab <- -log(-log(1 - 1 / Ttick_lab))

va_range <- range(Q)
ymax     <- max(UB) * 1.05

plot(
  y_emp, Q,
  xaxt = "n",
  xlab = "Return period T (years)",
  ylab = "Annual peak discharge (cfs)",
  main = paste("Gumbel frequency curve \u2014 USGS", mysite),
  pch  = 21, bg = "#8aad96", col = "#1e3a2f", cex = 1.1,
  ylim = c(0, ymax),
  xlim = c(min(ytick_fine), max(ytick_fine))
)
axis(1, at = ytick_lab, labels = as.character(Ttick_lab))

# Fitted line and confidence limits
lines(ytick_fine, Q_T, col = "#1a1a18", lwd = 2)
lines(ytick_fine, LB, col = "#2e7d94", lwd = 1.5, lty = 2)
lines(ytick_fine, UB, col = "#c0392b", lwd = 1.5, lty = 2)

# Grid and reference line
abline(v = ytick_fine, lty = 3, col = "lightgray")
abline(h = seq(0, ymax, length.out = 6), lty = 3, col = "lightgray")
abline(h = Qmax_ref, lty = 2, col = "#c0392b", lwd = 2)

legend("topleft",
       legend = c("Annual peaks", "Gumbel fit",
                  "95% CI", paste("Ref Q =", format(Qmax_ref, big.mark = ","))),
       pch  = c(21, NA, NA, NA),
       lty  = c(NA, 1,  2,  2),
       lwd  = c(NA, 2, 1.5, 2),
       col  = c("#1e3a2f", "#1a1a18", "#2e7d94", "#c0392b"),
       pt.bg = c("#8aad96", NA, NA, NA),
       bty = "n", cex = 0.82)

6.1 Estimated quantiles at standard return periods

T_standard <- c(2, 5, 10, 25, 50, 100)
y_std   <- -log(-log(1 - 1 / T_standard))
K_std   <- -(sqrt(6) / pi) * (0.5772 + log(log(T_standard / (T_standard - 1))))
Q_std   <- xbar + K_std * s
se_std  <- (s * sqrt(1 + 1.14 * K_std + 1.1 * K_std^2)) / sqrt(nQ)
lb_std  <- Q_std - qt(0.975, nQ - 1) * se_std
ub_std  <- Q_std + qt(0.975, nQ - 1) * se_std

tibble(
  T              = T_standard,
  Q_T            = round(Q_std),
  Lower_95       = round(lb_std),
  Upper_95       = round(ub_std),
  exceed_prob    = round(1 / T_standard, 4)
) %>%
  kbl(
    caption   = paste("Gumbel flood quantiles \u2014 USGS", mysite),
    col.names = c("Return period (yr)", "Q_T (cfs)",
                  "Lower 95% CI", "Upper 95% CI", "Annual exceedance prob.")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)
Gumbel flood quantiles — USGS 07061000
Return period (yr) Q_T (cfs) Lower 95% CI Upper 95% CI Annual exceedance prob.
2 109344 93994 124694 0.50
5 181897 156044 207750 0.20
10 229934 195015 264852 0.10
25 290628 243546 337709 0.04
50 335654 279317 391991 0.02
100 380348 314720 445976 0.01

7. Log-Pearson Type III

The Log-Pearson Type III (LP3) distribution is the standard method recommended by the US Water Resources Council (Bulletin 17C) for federal flood frequency studies. It fits a Pearson III distribution to the base-10 logarithms of annual peaks, capturing the positive skewness typical of flood series.

# Step 1: log-transform
logQ      <- log10(Q)

# Step 2: statistics of log-values
mean_logQ <- mean(logQ)
sd_logQ   <- sd(logQ)
skew_logQ <- skewness(logQ)

cat(sprintf("Log10(Q): mean = %.4f  SD = %.4f  skewness = %.4f\n",
            mean_logQ, sd_logQ, skew_logQ))
## Log10(Q): mean = 5.0043  SD = 0.2723  skewness = 0.0972
# Step 3: return periods and Pearson III frequency factor
T_lp3 <- c(2, 5, 10, 25, 50, 100)
z     <- qnorm(1 - 1 / T_lp3)          # standard normal variate

# Wilson-Hilferty approximation for Pearson III frequency factor K
g     <- skew_logQ
K_lp3 <- z + (z^2 - 1) * g/6 - (z^3 - 6*z) * g^2/36 +
         (z^2 - 1) * g^3/216 + z * g^4/648

# Step 4: estimated flood quantiles
logQ_est <- mean_logQ + K_lp3 * sd_logQ
Q_lp3    <- 10^logQ_est

tibble(
  T           = T_lp3,
  K           = round(K_lp3, 3),
  log10_Q_est = round(logQ_est, 4),
  Q_est_cfs   = round(Q_lp3)
) %>%
  kbl(
    caption   = paste("Log-Pearson III flood quantiles \u2014 USGS", mysite),
    col.names = c("Return period (yr)", "K factor",
                  "log\u2081\u2080(Q) estimate", "Q estimate (cfs)")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)
Log-Pearson III flood quantiles — USGS 07061000
Return period (yr) K factor log₁₀(Q) estimate Q estimate (cfs)
2 -0.016 4.9999 99978
5 0.838 5.2325 170822
10 1.293 5.3566 227276
25 1.786 5.4906 309427
50 2.107 5.5781 378507
100 2.398 5.6574 454378

📌 Gumbel vs. Log-Pearson III. Both methods produce similar estimates near the median (2–10 year return periods). They diverge at the extremes: LP3 accommodates positive skewness in the log domain more flexibly, which typically produces higher estimates for rare events (100-year floods) when the record has heavy-tailed behaviour. If the two methods disagree substantially at T = 100, inspect the Q-Q plots for both distributions before choosing one to report.

8. Stationarity Check — Two-Period Comparison

Splitting the record and comparing frequency curves is a visual diagnostic for non-stationarity. A rightward shift in the curve (higher discharge for the same return period) in the recent period suggests an intensification of flood magnitudes.

period1 <- annualpeak %>%
  filter(peak_dt >= min(peak_dt) & peak_dt < as.Date(splitdate))

period2 <- annualpeak %>%
  filter(peak_dt >= as.Date(splitdate))

cat(sprintf("Period 1: %d peaks  (%s to %s)\n",
            nrow(period1),
            format(min(period1$peak_dt), "%Y"),
            format(max(period1$peak_dt), "%Y")))
## Period 1: 71 peaks  (1904 to 1994)
cat(sprintf("Period 2: %d peaks  (%s to %s)\n",
            nrow(period2),
            format(min(period2$peak_dt), "%Y"),
            format(max(period2$peak_dt), "%Y")))
## Period 2: 24 peaks  (2001 to 2024)
# Helper: fit Gumbel and return a list of plotting ingredients
fit_gumbel <- function(df) {
  Q    <- df$peak_va
  n    <- length(Q)
  xbar <- mean(Q);  s <- sd(Q)
  r    <- n + 1 - rank(Q)
  T_e  <- (n + 1) / r
  y_e  <- -log(-log(1 - 1 / T_e))

  Ttick <- c(1.001,1.01,1.1,1.5,2,3,4,5,6,7,8,9,10,15,20,25,
             30,35,40,50,60,70,80,90,100)
  yt    <- -log(-log(1 - 1 / Ttick))
  K     <- -(sqrt(6)/pi) * (0.5772 + log(log(Ttick/(Ttick-1))))
  Q_T   <- xbar + K * s
  se    <- (s * sqrt(1 + 1.14*K + 1.1*K^2)) / sqrt(n)
  LB    <- Q_T - qt(0.975, n-1)*se
  UB    <- Q_T + qt(0.975, n-1)*se

  graphlab <- paste(range(format(df$peak_dt, "%Y")), collapse = " \u2013 ")
  list(Q=Q, n=n, y_emp=y_e, T_emp=T_e,
       ytick=yt, Q_T=Q_T, LB=LB, UB=UB, graphlab=graphlab)
}

g1 <- fit_gumbel(period1)
g2 <- fit_gumbel(period2)
ymax_all <- max(c(g1$UB, g2$UB)) * 1.05
xlims    <- c(min(g1$ytick), max(g1$ytick))

Ttick_lab <- c(1.01, 1.1, 1.5, 2, 5, 10, 20, 50, 100)
ytick_lab <- -log(-log(1 - 1 / Ttick_lab))

par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3.5, 1))

for (g in list(g1, g2)) {
  plot(
    g$y_emp, g$Q,
    xaxt = "n",
    xlab = "Return period T (years)",
    ylab = "Annual peak discharge (cfs)",
    main = paste("USGS", mysite, "\n", g$graphlab),
    pch  = 21, bg = "#8aad96", col = "#1e3a2f", cex = 1.0,
    ylim = c(0, ymax_all),
    xlim = xlims,
    cex.main = 0.92
  )
  axis(1, at = ytick_lab, labels = as.character(Ttick_lab), cex.axis = 0.8)
  lines(g$ytick, g$Q_T, col = "#1a1a18", lwd = 2)
  lines(g$ytick, g$LB,  col = "#2e7d94", lwd = 1.4, lty = 2)
  lines(g$ytick, g$UB,  col = "#c0392b", lwd = 1.4, lty = 2)
  abline(v = ytick_lab, lty = 3, col = "lightgray")
  abline(h = seq(0, ymax_all, length.out = 6), lty = 3, col = "lightgray")
  abline(h = Qmax_ref, lty = 2, col = "#c0392b", lwd = 2)
}

par(mfrow = c(1, 1))
# Compare T-year quantiles between periods
T_cmp  <- c(2, 5, 10, 25, 50, 100)
yt_cmp <- -log(-log(1 - 1 / T_cmp))

mk_Q <- function(g) {
  xbar <- mean(g$Q); s <- sd(g$Q)
  K    <- -(sqrt(6)/pi) * (0.5772 + log(log(T_cmp/(T_cmp-1))))
  round(xbar + K * s)
}

tibble(
  T      = T_cmp,
  Q_p1   = mk_Q(g1),
  Q_p2   = mk_Q(g2),
  change = paste0(ifelse(mk_Q(g2) >= mk_Q(g1), "+", ""),
                  round(100 * (mk_Q(g2) - mk_Q(g1)) / mk_Q(g1), 1), "%")
) %>%
  kbl(
    caption   = paste("Gumbel flood quantiles: Period 1 vs. Period 2 \u2014 USGS", mysite),
    col.names = c("Return period (yr)",
                  paste0("Q_T Period 1 (", g1$graphlab, ")"),
                  paste0("Q_T Period 2 (", g2$graphlab, ")"),
                  "Change (%)")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  column_spec(4, bold = TRUE,
              color = ifelse(
                mk_Q(g2) >= mk_Q(g1), "#c0392b", "#2e7d32"
              ))
Gumbel flood quantiles: Period 1 vs. Period 2 — USGS 07061000
Return period (yr) Q_T Period 1 (1904 – 1994) Q_T Period 2 (2001 – 2024) Change (%)
2 119620 81716 -31.7%
5 197828 122629 -38%
10 249608 149716 -40%
25 315033 183942 -41.6%
50 363568 209333 -42.4%
100 411746 234536 -43%

⚠️ Small-sample caution for split-period analysis. Splitting a 50-year record into two 25-year halves leaves each sub-period with very wide confidence intervals — especially at rare return periods (50, 100 year) that may not even be well-sampled within a 25-year window. Visual differences between the two curves are suggestive, not conclusive. Formal non-stationarity tests (e.g., Mann-Kendall trend test on the annual peak series, or a likelihood-ratio test on the Gumbel location parameter) are needed to make a defensible claim of change.

9. Reusable Helper Functions

The analysis above is packaged into three functions you can source into any project. They wrap the complete Gumbel pipeline and allow interactive classroom use.

getAnnualPeak.fun()

Downloads and cleans USGS annual peak data.

getAnnualPeak.fun <- function(site, startDate = "", endDate = "") {
  library(dataRetrieval)
  ap <- readNWISpeak(site, startDate, endDate)
  ap <- ap[!is.na(ap$peak_va) & !is.na(ap$peak_dt), ]
  ap
}

splitdata.fun()

Splits an annual peak data frame at a given date.

splitdata.fun <- function(annualpeak, splitdate) {
  p1 <- subset(annualpeak, peak_dt <  as.Date(splitdate))
  p2 <- subset(annualpeak, peak_dt >= as.Date(splitdate))
  list(period1 = p1, period2 = p2)
}

plot_floodreturn.fun()

Plots the Gumbel frequency curve for any annual peak data frame.

plot_floodreturn.fun <- function(annualpeak, Qmax_event = NULL) {

  Q        <- annualpeak$peak_va
  va_range <- range(Q)
  site_no  <- annualpeak$site_no[1]
  graphlab <- paste(range(format(annualpeak$peak_dt, "%Y")), collapse = " \u2013 ")

  n     <- length(Q)
  r     <- n + 1 - rank(Q)
  T_emp <- (n + 1) / r
  y_emp <- -log(-log(1 - 1 / T_emp))

  Ttick <- c(1.001, 1.01, 1.1, 1.5, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
             14, 15, 16, 17, 18, 19, 20, 25, 30, 35, 40, 45, 50, 60, 70,
             80, 90, 100)
  xtlab <- c(1.001, 1.01, 1.1, 1.5, 2, NA, NA, 5, NA, NA, NA, NA, 10, NA,
             NA, NA, NA, 15, NA, NA, NA, NA, 20, NA, 30, NA, NA, NA, 50,
             NA, NA, NA, NA, 100)
  ytick <- -log(-log(1 - 1 / Ttick))

  xbar  <- mean(Q);  s <- sd(Q)
  K     <- -(sqrt(6)/pi) * (0.5772 + log(log(Ttick/(Ttick-1))))
  Q_T   <- xbar + K * s
  se    <- (s * sqrt(1 + 1.14*K + 1.1*K^2)) / sqrt(n)
  LB    <- Q_T - qt(0.975, n - 1) * se
  UB    <- Q_T + qt(0.975, n - 1) * se
  ymax  <- if (!is.null(Qmax_event)) max(max(UB), Qmax_event) * 1.05
            else max(UB) * 1.05

  plot(
    y_emp, Q,
    ylab  = "Annual peak discharge (cfs)",
    xaxt  = "n",
    xlab  = "Return period T (years)",
    ylim  = c(0, ymax),
    xlim  = c(min(c(min(y_emp), min(ytick))), max(ytick)),
    pch   = 21, bg = "#8aad96", col = "#1e3a2f",
    main  = paste(site_no, graphlab),
    cex.main = 0.9
  )
  par(cex = 0.65)
  axis(1, at = ytick, labels = as.character(xtlab))
  par(cex = 1)

  lines(ytick, Q_T, col = "#1a1a18", lty = 1, lwd = 2)
  lines(ytick, LB,  col = "#2e7d94", lty = 2, lwd = 1.5)
  lines(ytick, UB,  col = "#c0392b", lty = 2, lwd = 1.5)

  abline(v = ytick, lty = 3, col = "lightgray")
  abline(h = seq(0, ymax, length.out = 6), lty = 3, col = "lightgray")

  if (!is.null(Qmax_event))
    abline(h = Qmax_event, lty = 2, col = "#c0392b", lwd = 2)
}

Quick demo of the functions

# Replicate the full-record plot using the helper function
plot_floodreturn.fun(annualpeak, Qmax_ref)

# Split and compare using the helper functions
periods <- splitdata.fun(annualpeak, splitdate)
par(mfrow = c(1, 2))
plot_floodreturn.fun(periods$period1, Qmax_ref)
plot_floodreturn.fun(periods$period2, Qmax_ref)

par(mfrow = c(1, 1))

10. Summary

💡 Key takeaways

💬 Discussion questions

  1. The April 2025 event at the Current River registered approximately 180,000 cfs. Based on your Gumbel frequency curve, what is the estimated return period of that event? How does this compare between the two sub-periods? What does the difference tell you?
  2. The Weibull formula cannot estimate return periods longer than n+1 years (where n is the record length). For a 75-year record, what is the longest return period you could estimate empirically — and why is extrapolating to T = 500 risky?
  3. The 95% confidence intervals are very wide for the 100-year quantile. Explain mathematically why rare quantiles have higher uncertainty than common quantiles, using the standard error formula in the Gumbel fit.
  4. Milly et al. (2008) argued that stationarity is dead. Describe two concrete mechanisms — one climatic, one anthropogenic — that could make the flood frequency distribution at the Current River non-stationary. How would each mechanism shift the fitted frequency curve?
  5. The Log-Pearson III and Gumbel methods will give similar results near the median but may diverge at T = 100. Run both for the full record at your station and compare. Which gives a higher 100-year estimate, and can you explain why based on the skewness of the log-transformed series?

References

Bulletin 17C. 2019. Guidelines for Determining Flood Flow Frequency. US Interagency Committee on Water Data, Hydrology Subcommittee. USGS Techniques and Methods 4-B5.

Milly, P. C. D., et al. 2008. Stationarity is dead: Whither water management? Science 319:573–574.

Merz, B., et al. 2021. Causes, impacts and patterns of disastrous river floods. Nature Reviews Earth and Environment 2:592–609.

USGS. 2019. National Water Information System. https://waterdata.usgs.gov/nwis

EA030 RTricks  |  Prof. Marc Los Huertos  |  Pomona College  |  Activity 11  |  Back to index