Recession Analysis

HydroRaVENS includes tools for analysing the nonlinear storage–discharge relationship directly from observed streamflow records, before running or calibrating the model.

Theory

During a baseflow recession (no rainfall or snowmelt), water drains from storage at a rate set by the storage–discharge relationship. If discharge is a power law of storage depth,

\[Q = c \, H^b,\]

then substituting \(\mathrm{d}H/\mathrm{d}t = -Q\) gives

\[-\frac{\mathrm{d}Q}{\mathrm{d}t} = K \, Q^n, \qquad n = \frac{2b - 1}{b},\]

where \(b\) is the HydroRaVENS recession_exponent and \(n\) is the slope on a log–log plot of \(-\mathrm{d}Q/\mathrm{d}t\) versus \(Q\) — the Brutsaert & Nieber (1977) recession plot.

Inverting:

\[b = \frac{1}{2 - n}\]

The two exponents have the following special values:

B–N slope n	HydroRaVENS b	Physical interpretation
1	1	Linear reservoir; exponential recession.
3/2	2	Long-time Boussinesq solution for horizontal, unconfined aquifer (Brutsaert & Nieber 1977).
→ 2	→ ∞	Extreme nonlinearity; storage empties nearly instantaneously.

The conversion \(b = 1/(2-n)\) is valid only for \(n < 2\). The Boussinesq short-time solution gives \(n = 3\), which is outside this framework; in practice, early fast-response recessions (surface runoff, tile drains) produce steep tails on the B–N cloud that should not be used to set the baseflow recession exponent.

Brutsaert–Nieber Analysis

class hydroravens.BrutsaertNieber(Q, dt=1.0, min_recession_days=3)

Brutsaert & Nieber (1977) recession analysis.

Extracts all recession segments from an observed specific-discharge record, computes –dQ/dt vs Q for each consecutive pair within a segment, and fits a power law on the resulting log–log cloud:

-dQ/dt = a * Q^n

The fitted slope n can be converted to a HydroRaVENS recession_exponent b via to_reservoir_exponent().

Parameters:

• Q (array-like) – Observed specific discharge time series [mm/day, or any consistent units]. Must be non-negative; zero values are skipped.

• dt (float, optional) – Timestep length [days]. Default 1.0.

• min_recession_days (int, optional) – Minimum number of consecutive declining days required for a segment to be included in the fit. Default 3.

n_

Fitted B–N slope (exponent in –dQ/dt = a·Q^n). Set by fit().

Type:

float

a_

Fitted coefficient. Set by fit().

Type:

float

r2_

Coefficient of determination of the power-law (linear log–log) fit. Set by fit(). Values below ~0.4 indicate a poor fit; a warning is issued automatically.

Type:

float

r2_quadratic_

R² of a degree-2 polynomial fit in log–log space. Set by fit(). When substantially larger than r2_, the recession cloud is systematically curved, suggesting the power-law storage–discharge assumption does not hold; a warning is issued automatically.

Type:

float

Examples

>>> import numpy as np
>>> from hydroravens import BrutsaertNieber
>>> Q = np.array([10.0, 8.0, 6.5, 5.2, 4.0, 3.1, 2.4, 1.8, 1.3])
>>> bn = BrutsaertNieber(Q).fit()
>>> print(f"n = {bn.n_:.3f},  b_HR = {bn.to_reservoir_exponent():.3f}")

fit()

Extract recession pairs, fit the log–log power law, and assess fit quality.

Performs ordinary least-squares regression of log(–dQ/dt) on log(Q). The resulting slope is stored in n_ and the coefficient in a_.

Two goodness-of-fit diagnostics are computed and stored:

• r2_ — R² of the linear (power-law) log–log fit.

• r2_quadratic_ — R² of a degree-2 polynomial log–log fit.

A UserWarning is issued if the recession cloud shows systematic curvature (quadratic fit substantially better than linear, ΔR² > 0.05), or if the power-law fit is poor overall (R² < 0.4). Curvature indicates that the power-law storage–discharge assumption may not hold; following Kirchner (2009), a non-parametric characterisation of the catchment sensitivity function should be considered before committing to a reservoir functional form.

Returns:

Allows method chaining (e.g. BrutsaertNieber(Q).fit()).

Return type:

self

Raises:

ValueError – If fewer than three recession pairs are found (cannot fit).

to_reservoir_exponent()

Convert the fitted B–N slope n to a HydroRaVENS recession_exponent b.

The relationship derives from substituting Q ∝ H^b and dH/dt = –Q into the B–N form –dQ/dt = a·Q^n:

n = (2b - 1) / b   →   b = 1 / (2 - n)

This conversion is valid only for n < 2. For n ≥ 2, the equivalent power-law reservoir exponent is infinite (or ill-defined), and a UserWarning is issued.

Returns:

HydroRaVENS recession_exponent b. Returns numpy.inf when n ≥ 2.

Return type:

float

Raises:

RuntimeError – If fit() has not been called.

summary()

Print a short summary of the fitted parameters.

Raises:

RuntimeError – If fit() has not been called.

plot(ax=None, show_fit=True, label=None, **scatter_kwargs)

Plot –dQ/dt vs Q on a log–log scale.

Parameters:

• ax (matplotlib.axes.Axes, optional) – Axes to draw on. Created if not supplied.

• show_fit (bool, optional) – Overlay the fitted power-law line. Default True.

• label (str, optional) – Label for the scatter points (appears in legend when show_fit=True).

• **scatter_kwargs – Passed to ax.scatter.

Return type:

matplotlib.axes.Axes

Raises:

RuntimeError – If fit() has not been called.

Using Results in HydroRaVENS

The B–N analysis produces a single catchment-integrated exponent b. In a multi-reservoir HydroRaVENS model, the reservoirs have different physical roles and the exponent should be assigned accordingly:

Reservoir	Suggested b	Rationale
Soil (fastest)	B–N estimate	The full recession cloud is dominated by the soil/fast response. Use as a calibration starting point; typical calibrated values are b ≈ 3–4 for agricultural catchments.
Intermediate	B–N lower-envelope estimate (fix, do not calibrate)	The lower envelope of the B–N cloud represents the slowest observable recession pathway — typically the shallow Quaternary zone (outwash, fractured regolith; MRT days to weeks). The truly deep reservoir drains too slowly and at too low a flux to register in the cloud, so B–N is capturing the intermediate, not the deep. The long-time Boussinesq value (\(b = 2.0\)) is a reasonable lower bound; basin-specific fits typically yield \(b \approx 2.0\)–2.5. Fix rather than calibrate (see equifinality note in Model Description).
Deep (slowest)	1.0 (fixed)	Mean residence time of decades to centuries; contributes too little flux to constrain b from streamflow. The linear approximation (Darcy flow in a confined or semi-confined aquifer) is physically appropriate and avoids adding a free parameter.

The B–N slope n from the lower envelope of the cloud corresponds to long-duration recessions dominated by the intermediate subsurface zone and is the most physically meaningful region for setting the intermediate-reservoir exponent. The upper scatter reflects short, event-driven recessions (tile drains, surface runoff) and should not be used to set the intermediate exponent.

For a complete prior-estimation workflow that combines the B–N analysis with timescale estimation, see suggest_priors() and the Tutorial: From Data to Calibrated Model.

Workflow

A typical workflow fits the recession curve and uses the result as a prior for model calibration:

import pandas as pd
from hydroravens import BrutsaertNieber

df = pd.read_csv('streamflow.csv', parse_dates=['Date'])
Q  = df['Specific Discharge [mm/day]'].values

bn = BrutsaertNieber(Q, min_recession_days=3).fit()
bn.summary()
bn.plot()

b_prior = bn.to_reservoir_exponent()
print(f"Suggested recession_exponent prior: {b_prior:.2f}")

The output might look like:

Brutsaert–Nieber recession analysis
  Recession pairs used : 412
  Fitted slope  n      : 1.5831
  Fitted coeff  a      : 0.0147
  HydroRaVENS   b      : 2.367
  Reference (long-time Boussinesq): n = 1.5, b = 2.0

The lower envelope of the B–N cloud corresponds to long-duration recessions dominated by the intermediate subsurface zone (shallow Quaternary units, fractured regolith; MRT days to weeks) and is the most physically meaningful region for setting the intermediate-reservoir exponent. The truly deep reservoir (MRT decades–centuries) produces discharge too small to dominate any part of the cloud and should be fixed at \(b = 1\) independently. The upper scatter reflects short, event-driven recessions (tile drains, surface runoff) that are better captured by a calibrated soil-zone exponent.

Goodness of fit and the power-law assumption

fit() automatically assesses whether the recession cloud is well described by a power law and issues a UserWarning when it is not.

Two diagnostics are computed and stored as attributes:

• R² of the power-law fit (r2_) — coefficient of determination of the linear log–log regression. The B–N cloud is inherently noisy (all types of recession contribute), so R² values of 0.4–0.7 are typical for real catchments. R² < 0.4 indicates a poor overall fit; inspect the plot.

• Curvature test (r2_quadratic_) — R² of a degree-2 polynomial fit in log–log space. If the quadratic fit is substantially better than the linear fit (ΔR² > 0.05), the recession cloud is systematically curved. A straight line in log–log space is the signature of a power-law storage–discharge relationship; systematic curvature means this assumption does not hold.

What to do when the power-law assumption fails:

A curved recession cloud in log–log space signals that the catchment sensitivity function g(Q) — as defined by Kirchner (2009) — is not a power law. This can arise from:

• Multiple reservoirs draining simultaneously, each with different exponents, producing a composite curve. Consider whether a multi-reservoir model with different fixed exponents captures the cloud’s shape better than a single-exponent fit.

• A genuine threshold or non-power-law drainage process (e.g. a perched water table that drains rapidly when full and slowly when nearly empty).

• Data quality issues (precipitation misattribution, rating curve errors at low flow).

In the first case, the lower envelope still provides a useful prior for the intermediate reservoir exponent. In the second case, a different reservoir functional form may be needed — one that HydroRaVENS does not currently implement, following Kirchner’s (2009) more general framework. In the third case, inspect and correct the input data before re-fitting.

Interpretation guide

• n ≈ 1.5, b ≈ 2 — recession consistent with the Boussinesq long-time solution for a horizontal unconfined aquifer. A reasonable fixed prior for the intermediate reservoir (shallow Quaternary units, fractured regolith). The deep reservoir should be fixed at \(b = 1\) regardless of the B–N result.

• n ≈ 1.3–1.6, b ≈ 1.5–2.5 — typical mixed catchment range.

• n > 1.7, b > 3 — rapid nonlinear drainage, often reflecting near-surface soil processes. Calibrate rather than fix.

• n ≥ 2 — unphysical in this framework; refit using only the lower envelope or with a longer min_recession_days.

References

Brutsaert, W. and Nieber, J. L. (1977). Regionalized drought flow hydrographs from a mature glaciated plateau. Water Resources Research, 13(3), 637–643. https://doi.org/10.1029/WR013i003p00637

Kirchner, J. W. (2009). Catchments as simple dynamical systems: Catchment characterization, rainfall-runoff modeling, and doing hydrology backward. Water Resources Research, 45, W02429. https://doi.org/10.1029/2008WR006912