Flood Frequency and Streamflow Statistics

Where a stream gauge has been recording for decades, the flood of a stated recurrence interval can be estimated from what the river has actually done, with no rainfall model, no loss method and no unit hydrograph in the chain. As part of the rainfall-runoff modeling and hydrologic simulation workflow, frequency analysis is both an alternative to the modelling chain and the thing that chain must be checked against — a modelled 100-year peak that disagrees with the gauged frequency curve at the same site is a finding about the model, not about the river.

It is also the part of hydrology where the honest answer is most often “wider than you would like”. A frequency curve fitted to twenty years of record and extrapolated to a hundred-year return period carries a confidence interval that spans a factor of two, and reporting the central estimate without it is the most common way these numbers mislead.

Prerequisites and Environment Setup

conda create -n floodfreq python=3.11
conda activate floodfreq
conda install -c conda-forge pandas=2.2 numpy scipy=1.12 matplotlib
pip install dataretrieval==1.0.9

The qualification codes matter more than they look. A peak flagged as affected by regulation is not a natural flood and cannot be pooled with natural ones; a peak flagged as a historic estimate carries different weight and, strictly, calls for a historically-weighted fit. Silently dropping the codes on retrieval is a common and consequential shortcut.

Mechanics: What a Frequency Curve Claims

A flood frequency curve says: if the future resembles the record, a flow of this size will be exceeded in this fraction of years. Three assumptions are buried in that sentence, and each fails in identifiable ways.

Why Log-Pearson III, and what the skew does

Annual peaks are right-skewed: most years cluster low and a few sit far above them. Taking logarithms brings the distribution close to symmetric, and fitting a Pearson Type III — a gamma family with a shift — to those logarithms adds one parameter, the skew, that controls how heavy the upper tail is.

That single parameter dominates the extrapolation. The mean and standard deviation are estimated well from a short record; the skew is not, because it depends on the third moment and is therefore extremely sensitive to the largest one or two peaks. This is why the standard practice weights the station skew against a regional value.

A station skew estimated from 25 years has a standard error around 0.4, which is the entire range of that table. Weighting against a regional skew with a smaller mean square error is not a formality; it is what stops one exceptional flood from setting the design value.

Step-by-Step Workflow

1. Retrieve the annual peak series with qualification codes. Keep the water-year convention consistent — mixing calendar and water years splits single events across two years.
2. Screen the record. Drop regulated years or fit them separately. Test for trend and for a change point. Test lag-one autocorrelation.
3. Fit Log-Pearson III to the base-ten logarithms: compute the mean, standard deviation and skew of the logs.
4. Weight the skew against the regional value, using the inverse of each estimate’s mean square error as its weight.
5. Compute quantiles for the return periods required, using the frequency factor for the weighted skew.
6. Attach confidence limits derived from the record length, and report them alongside the estimates rather than in an appendix.
7. Sanity-check against neighbouring gauges by unit discharge, and against the delineated drainage area.

Production-Ready Code

The function below fits Log-Pearson III with a weighted skew and returns quantiles with confidence limits.

import logging

import numpy as np
import pandas as pd
from scipy import stats

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")


def fit_log_pearson3(
    annual_peaks: pd.Series,
    regional_skew: float | None = None,
    regional_skew_mse: float = 0.302,
    return_periods: tuple[int, ...] = (2, 5, 10, 25, 50, 100, 200, 500),
    confidence: float = 0.90,
) -> pd.DataFrame:
    """
    Fit a Log-Pearson Type III flood frequency curve to an annual peak series.

    Parameters
    ----------
    annual_peaks      : One peak per water year, indexed by year. Zero and
                        negative values are dropped with a warning.
    regional_skew     : Published regional skew. When given, it is weighted
                        against the station skew by inverse mean square error.
    regional_skew_mse : Mean square error of the regional skew value.
    return_periods    : Return periods to report, in years.
    confidence        : Two-sided confidence level for the reported limits.

    Returns
    -------
    DataFrame indexed by return period with columns: exceedance_prob,
    estimate, lower, upper.
    """
    peaks = pd.Series(annual_peaks).dropna()
    n_all = len(peaks)
    peaks = peaks[peaks > 0]
    if len(peaks) < n_all:
        log.warning("Dropped %d non-positive peak(s) — a zero annual peak needs "
                    "conditional-probability treatment, not deletion",
                    n_all - len(peaks))

    n = len(peaks)
    if n < 10:
        raise ValueError(f"only {n} usable years — too short for a frequency fit")
    if n < 30:
        log.warning("Record is %d years; estimates beyond the %d-year return "
                    "period are extrapolation with wide limits", n, 2 * n)

    logs = np.log10(peaks.to_numpy(dtype=float))
    mean_log = float(logs.mean())
    std_log = float(logs.std(ddof=1))
    station_skew = float(stats.skew(logs, bias=False))

    # --- Weight the station skew against the regional value. The station skew
    # is the least reliable of the three moments, and its MSE grows sharply as
    # the record shortens, so this weighting is doing real work. ---
    skew = station_skew
    if regional_skew is not None:
        # Bulletin 17-style approximation of the station skew MSE.
        a = -0.33 + 0.08 * abs(station_skew) if abs(station_skew) <= 0.9 else \
            -0.52 + 0.30 * abs(station_skew)
        b = 0.94 - 0.26 * abs(station_skew) if abs(station_skew) <= 1.5 else 0.55
        station_mse = 10 ** (a - b * np.log10(n / 10.0))
        w = regional_skew_mse / (regional_skew_mse + station_mse)
        skew = w * station_skew + (1 - w) * regional_skew
        log.info("Station skew %.3f (MSE %.3f), regional %.3f (MSE %.3f) → "
                 "weighted %.3f (station weight %.2f)",
                 station_skew, station_mse, regional_skew, regional_skew_mse, skew, w)
    else:
        log.info("Station skew %.3f used unweighted — supply a regional skew "
                 "for a defensible design estimate", station_skew)

    z_conf = stats.norm.ppf(0.5 + confidence / 2.0)
    rows = []
    for T in return_periods:
        p_exceed = 1.0 / T
        # Frequency factor for Pearson III via the Wilson-Hilferty transform.
        z = stats.norm.ppf(1.0 - p_exceed)
        k = skew / 6.0
        kt = (2.0 / skew) * (((z - k) * k + 1.0) ** 3 - 1.0) if abs(skew) > 1e-6 else z

        log_q = mean_log + kt * std_log
        estimate = 10 ** log_q

        # Standard error of the fitted quantile grows with Kt and shrinks with n.
        se_log = std_log * np.sqrt((1.0 + 0.5 * kt ** 2) / n)
        lower = 10 ** (log_q - z_conf * se_log)
        upper = 10 ** (log_q + z_conf * se_log)

        rows.append({
            "return_period": T,
            "exceedance_prob": p_exceed,
            "estimate": estimate,
            "lower": lower,
            "upper": upper,
            "band_ratio": upper / lower,
        })
        log.info("T=%4d yr: %8.1f  [%.1f, %.1f]  band ×%.2f",
                 T, estimate, lower, upper, upper / lower)

    out = pd.DataFrame(rows).set_index("return_period")
    log.info("Fit from %d years: mean(log)=%.4f sd(log)=%.4f skew=%.3f",
             n, mean_log, std_log, skew)
    return out


def screen_record(annual_peaks: pd.Series) -> dict:
    """Run the three assumption tests and report, without deciding for you."""
    s = pd.Series(annual_peaks).dropna()
    years = np.asarray(s.index, dtype=float)
    vals = s.to_numpy(dtype=float)

    tau, p_trend = stats.kendalltau(years, vals)
    lag1 = float(pd.Series(vals).autocorr(lag=1))

    # A crude change-point scan: the split that maximises the difference in means.
    best_split, best_stat = None, 0.0
    for i in range(5, len(vals) - 5):
        t_stat = abs(vals[:i].mean() - vals[i:].mean()) / (vals.std(ddof=1) + 1e-9)
        if t_stat > best_stat:
            best_stat, best_split = t_stat, int(s.index[i])

    result = {
        "n_years": len(s),
        "kendall_tau": float(tau),
        "trend_p_value": float(p_trend),
        "lag1_autocorrelation": lag1,
        "change_point_year": best_split,
        "change_point_strength": float(best_stat),
    }
    if p_trend < 0.05:
        log.warning("Significant trend (tau=%.3f, p=%.4f) — stationarity is "
                    "questionable; consider a non-stationary fit or a split record",
                    tau, p_trend)
    if abs(lag1) > 0.3:
        log.warning("Lag-1 autocorrelation %.2f — peaks are not independent; "
                    "confidence limits will be optimistic", lag1)
    return result


# --- Example usage ---
# import dataretrieval.nwis as nwis
# peaks_df, _ = nwis.get_discharge_peaks(sites="03339000")
# peaks = peaks_df.set_index(peaks_df["peak_dt"].dt.year)["peak_va"]
# print(screen_record(peaks))
# print(fit_log_pearson3(peaks, regional_skew=-0.05).round(1))

Validation Protocol

• Plotting positions against the fit. Plot the ranked observations at their Weibull plotting positions on the same axes as the fitted curve. Systematic departure in the upper tail means the distribution choice is wrong for this record, not that the data is odd.
• Unit discharge against neighbours. Divide the 100-year estimate by drainage area and compare against nearby gauges on comparable basins. An outlier by more than a factor of two needs an explanation before the number is used.
• Leave-one-out sensitivity. Refit with the largest peak removed. If the 100-year estimate moves by more than about 20 %, the result is one flood’s opinion and the confidence band should be reported prominently.
• Area consistency. The published drainage area for the gauge should match the delineated catchment. A mismatch means the frequency curve and the modelled basin are not the same object — see snapping stream gauge locations to NHD flowlines.
• Cross-check against the modelled peak. A rainfall-runoff model driven by a design storm of the same return period should land inside the frequency curve’s confidence band. Outside it, the model needs attention.

Low-Flow Statistics

The same record answers questions at the other end of the distribution, and the mechanics differ enough to be worth stating.

7Q10 — the lowest seven-day mean flow with a ten-year recurrence — is computed from a rolling seven-day mean of daily flows, taking one annual minimum per climatic year, then fitting the lower tail. Three details separate a correct calculation from a plausible one:

• The climatic year for low flow runs April to March in most of the northern hemisphere, so a single low-flow season is not split across two years.
• Zero-flow years cannot be log-transformed. They need a conditional-probability treatment: fit the non-zero years and adjust the exceedance probabilities by the fraction of zero years.
• The fitted distribution is usually log-Pearson III again, but fitted to the minima, where the skew has the opposite sign and the tail of interest is the lower one.

Common Failure Modes and Optimization

• Extrapolating far past the record. A 20-year record does not contain information about a 500-year flood. The arithmetic produces a number; the number is not an estimate.
• Unweighted station skew. One exceptional flood can move the 100-year estimate by 30 %. Weighting against a regional skew is the standard defence and takes one extra argument.
• Mixed populations. Snowmelt and rainfall peaks pooled into one series produce a curve that fits neither and shows a visible kink in the plotting positions. Fit them separately and combine the exceedance probabilities.
• Regulated years included. A reservoir upstream truncates the upper tail after its construction date. Split the record at the construction year, and use only the natural period unless you are explicitly estimating regulated flows.
• Water year used for low flow. Splits a single drought across two years and understates the severity of both. Use the climatic year.
• Zero-flow years deleted. Deleting them biases the low-flow statistic upward, sometimes severely on ephemeral streams. Conditional probability is the correct treatment.
• Confidence limits omitted. A single number implies a precision that the fit does not have. Report the band; if the band is embarrassing, that is information.

When to Use This vs. Alternatives

Use frequency analysis wherever a long, natural, homogeneous record exists at or very near the site of interest. It is the most direct evidence available and involves no rainfall model at all.

Use regional regression at an ungauged site. Published regional equations relate quantiles to drainage area and a few climatic variables and are calibrated on many gauges at once. They are less precise at any one site but do not require a record.

Use a rainfall-runoff model when the catchment has changed, or is about to. A frequency curve describes the basin the record came from; if that basin is being urbanised or dammed, the record no longer describes it, and only a model can represent the change. The precipitation forcing and design storms page covers the input side of that chain.

Use both, and compare. Where a gauge exists, fitting the curve and running the model is not duplicated effort: agreement is evidence, and disagreement localises a problem in one of them.

Frequently Asked Questions

How long a gauge record is needed for a 100-year flood estimate?

There is no record length that makes a 100-year estimate precise, but the usual guidance is that a record should be at least a tenth of the target return period before the estimate is meaningful, and comfortably longer before it is stable. With 20 years the 100-year confidence interval typically spans a factor of two; with 60 years it narrows to roughly ±25 %. Extrapolating a 15-year record to a 500-year flood is arithmetic, not estimation.

Why is flood frequency fitted in log space?

Annual flood peaks are strongly right-skewed and strictly positive, and their logarithms are much closer to symmetric. Fitting a Pearson Type III distribution to the logarithms handles that skew with a single extra parameter, keeps every fitted quantile positive, and matches the convention that regional skew maps and published station statistics are expressed in.

What is 7Q10 and how is it different from a flood quantile?

7Q10 is the lowest seven-day average flow expected once in ten years — a low-flow statistic used for water-quality permitting and instream flow rules. It is computed from annual minima of a seven-day rolling mean rather than from annual maxima, and it is fitted in the lower tail, where the distribution behaves differently and zero-flow years need explicit handling. A flood quantile and a low-flow statistic come from the same record but from opposite ends of it.

Related Topics

• Rainfall-Runoff Modeling & Hydrologic Simulation — the parent section, and the modelling chain this analysis is the alternative to
• Precipitation Forcing and Design Storms — the design-storm route to the same quantity, for sites with no gauge
• Model Calibration & Objective Functions — where the gauge record is used to fit a model rather than a distribution
• Snapping Stream Gauge Locations to NHD Flowlines — confirming that the gauge and the delineated basin are the same catchment
• Outlet Point Mapping & Validation — the delineation side of that check