Unit Hydrograph Methods in Python

The unit hydrograph is the single most reused abstraction in event-based rainfall-runoff work: it collapses the entire routing behaviour of a basin into one transfer function that you convolve against excess rainfall to predict a storm hydrograph. Within the Rainfall-Runoff Modeling & Hydrologic Simulation domain, unit hydrograph methods sit between the loss calculation that produces excess rainfall and the channel routing that moves the resulting flood wave downstream. This guide treats the unit hydrograph as a linear, time-invariant (LTI) system, derives the discrete convolution equation, walks through the synthetic UH families you will actually reach for, and gives one production-grade convolution function in numpy. The excess-rainfall input this whole method consumes typically comes from a loss model such as the curve number runoff estimation approach, so the two topics are designed to be read together.

Two child guides drill into the parts that generate the transfer function itself: deriving SCS dimensionless unit hydrographs in Python for the single-parameter NRCS method, and Snyder vs Clark unit hydrograph: when to use each for the two dominant synthetic alternatives. The material below is the shared theory and the Python plumbing that all of them feed into.

The Unit Hydrograph as a Linear, Time-Invariant System

A unit hydrograph is the direct-runoff hydrograph (DRH) produced by exactly one unit of excess rainfall — 1 cm or 1 inch — generated uniformly over the whole basin at a constant rate during a specified duration D. Every practical use of the UH rests on three assumptions that together make the basin behave as an LTI system:

• Proportionality (linearity). If 1 unit of excess produces ordinates U(t), then 2 units produce 2·U(t). The response scales linearly with input depth.
• Superposition. The DRH from several excess-rainfall pulses is the sum of the individual responses, each lagged to its own pulse. This is what makes convolution valid.
• Time invariance. The UH is the same regardless of when the storm occurs. A pulse today and an identical pulse next week produce identically shaped responses, merely offset in time.

These assumptions are approximations. Real basins are mildly nonlinear — larger storms travel faster and peak proportionally higher — and antecedent moisture shifts the response. The UH is therefore an event-scale engineering model, not a physical law, but within a design-storm workflow it is accurate enough and computationally trivial. The two properties that define the UH numerically are its duration D (the excess pulse length that produced it), its lag (the delay from the centroid of excess rainfall to the peak), and its time to peak Tp (from the start of excess to the peak ordinate). Getting D consistent with your rainfall time step is the single most common source of error, and it is treated explicitly below.

Prerequisites & Environment Setup

The convolution mathematics needs almost nothing beyond numpy. The wider workflow — reading design storms, tabulating ordinates, exporting results — uses pandas, and scipy provides an alternate convolution primitive plus interpolation for scaling dimensionless tables.

conda create -n hydro-uh python=3.10 numpy scipy pandas -c conda-forge
conda activate hydro-uh

The inputs you must have in hand before writing any convolution code are an excess-rainfall hyetograph (depth per time step, already net of infiltration and other losses) and a unit hydrograph whose duration matches that time step. The excess hyetograph is the output of a loss method; if you are starting from gross rainfall, run it through curve number runoff estimation first so that the depths you convolve are genuinely excess and not gross precipitation.

Convolution Mechanics

Superposition of scaled, lagged unit hydrographs is exactly a discrete convolution. If the excess hyetograph has M pulses of depth P_m (m = 1…M) and the UH has J ordinates U_j (j = 1…J), the direct-runoff ordinate at time step n is:

Q_n = Σ  P_m · U_(n - m + 1)      for m from max(1, n-J+1) to min(n, M)

The resulting DRH has N = M + J − 1 ordinates. Each excess pulse P_m contributes a full copy of the UH, scaled by P_m and shifted so its leading edge starts at pulse m. The diagram below shows that decomposition: three excess pulses each scale and shift the same UH, and the vertical sum of those three responses is the composite direct-runoff hydrograph.

The duration/lag/time-to-peak relationships tie the shape to physical basin properties. For most synthetic families the time to peak follows Tp = D/2 + Tlag, where Tlag is the basin lag (centroid of excess to peak) and D is the UH duration. Because Tlag is a basin constant, changing D shifts Tp and changes the peak ordinate: a shorter, more intense unit pulse produces a higher, earlier peak. That coupling is precisely why the UH duration must match your rainfall time step, and why the S-curve exists to convert between durations.

Synthetic Unit Hydrograph Families

When a basin is gauged you can derive the UH directly by deconvolving observed storm hydrographs. Most basins are ungauged, so you build a synthetic UH from measurable basin characteristics. Four families dominate:

• SCS (NRCS) dimensionless UH. A single fixed dimensionless curve of t/Tp versus q/qp is scaled to the basin through Tp and a peak discharge qp = C·A/Tp. It needs only the basin lag (or time of concentration) and area, which makes it the default for ungauged design work. Full derivation is in deriving SCS dimensionless unit hydrographs in Python.
• Snyder synthetic UH. Basin lag tp = Ct·(L·Lca)^0.3 and peak qp = C·Cp·A/tp are driven by two regional coefficients Ct and Cp calibrated from gauged basins in the same physiographic region. It defines the peak, lag, and time base; the intervening shape is filled in with empirical width relations.
• Clark time-area UH. Routes rainfall through a time-area histogram (translation along isochrones) then through a single linear reservoir of storage coefficient R. Its two parameters, the time of concentration Tc and R, give a physically distributed response but require a GIS time-area curve.
• Nash cascade. Models the basin as n identical linear reservoirs in series, giving an instantaneous UH that is a gamma distribution: u(t) = (1/(k·Γ(n)))·(t/k)^(n−1)·e^(−t/k), with shape parameter n and storage constant k. It is compact and continuous, useful when you want a smooth analytic IUH to fit.

The choice between Snyder and Clark is a recurring decision that hinges on data availability and basin type; it has its own dedicated comparison in Snyder vs Clark unit hydrograph: when to use each. Whichever family you pick, the output is the same object — an array of UH ordinates at your time step — and everything downstream is convolution.

The S-Curve and Duration Conversion

The S-curve (or S-hydrograph) is the direct-runoff response to a continuous excess input applied at a constant rate of one unit per duration D, sustained indefinitely. You build it by summing copies of the D-hour UH, each lagged by D:

S(t) = U(t) + U(t − D) + U(t − 2D) + …

The S-curve climbs to an equilibrium plateau equal to the constant supply rate and is the tool for changing a UH’s duration. To convert a D-hour UH into a UH of a different duration D′, offset the S-curve by D′, subtract, and rescale by the ratio of durations:

U'(t) = (D / D') · [ S(t) − S(t − D') ]

This is essential in practice because published synthetic UHs and design hyetographs rarely share the same native time step. If your excess hyetograph is tabulated at 15-minute intervals but your UH was derived for a 1-hour duration, you must S-curve the UH down to 15 minutes before convolving. Skipping this step is the classic duration-mismatch error described under failure modes. The S-curve also exposes numerical noise: an S-curve that oscillates around its plateau instead of settling flat means the underlying UH ordinates are inconsistent with the duration, and any duration conversion built on it will inherit that ringing.

Step-by-Step Workflow

Step 1 — Compute the excess-rainfall hyetograph

Subtract losses from the gross design storm to leave excess rainfall at a fixed time step Δt. The SCS curve number method is the usual loss model; it returns cumulative excess depth, which you difference into per-step pulses. Every downstream ordinate inherits this time step, so pick it once and hold it constant.

import logging
import numpy as np

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

def incremental_excess(cumulative_excess_mm: np.ndarray) -> np.ndarray:
    """Difference a cumulative excess-depth curve into per-step excess pulses (mm)."""
    pulses = np.diff(cumulative_excess_mm, prepend=0.0)
    pulses = np.clip(pulses, 0.0, None)  # guard against tiny negative rounding
    logging.info("Excess pulses: %d steps, total depth %.2f mm", pulses.size, pulses.sum())
    return pulses

Step 2 — Build the unit hydrograph

Construct the UH ordinates for your chosen family, scaled so that the integral of the UH equals one unit of excess depth over the basin area. For the SCS method that means locating Tp and qp and sampling the dimensionless table; the numerical details live in the SCS dimensionless UH guide.

Step 3 — Match duration to the time step

Confirm the UH duration D equals Δt. If not, convert it through the S-curve as shown above. Convolving a 1-hour UH against 15-minute excess pulses silently multiplies your runoff volume by four.

Step 4 — Convolve

Apply the discrete convolution equation with numpy (Step covered in full by the production function below).

Step 5 — Add baseflow and export

Direct runoff is only the storm response. Add a baseflow component — a constant, a linear recession, or a master recession curve — to obtain total streamflow, then write the time-indexed ordinates for downstream routing or comparison against a gauge.

Production-Ready Code

The function below convolves an excess hyetograph with a unit hydrograph and returns a time-indexed direct-runoff hydrograph. It enforces the volume-conservation contract, checks duration consistency, guards against negative ordinates, and adds baseflow. This is the reusable core that every synthetic UH family feeds into.

import logging
import numpy as np
import pandas as pd

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")


def convolve_unit_hydrograph(
    excess_pulses_mm: np.ndarray,
    uh_ordinates: np.ndarray,
    dt_hours: float,
    area_km2: float,
    baseflow_cms: float = 0.0,
) -> pd.DataFrame:
    """
    Convolve an excess-rainfall hyetograph with a unit hydrograph.

    Parameters
    ----------
    excess_pulses_mm : per-step excess-rainfall depths (mm), one value per time step
    uh_ordinates     : unit hydrograph ordinates (m3/s per mm of excess), UH duration == dt_hours
    dt_hours         : computation time step in hours (must equal the UH duration)
    area_km2         : basin drainage area in square kilometres
    baseflow_cms     : constant baseflow added to direct runoff (m3/s)

    Returns
    -------
    DataFrame indexed by time (hours) with direct_runoff_cms and total_flow_cms columns.
    """
    excess = np.asarray(excess_pulses_mm, dtype=float)
    uh = np.asarray(uh_ordinates, dtype=float)

    # --- Guard: negative ordinates are physically impossible ---
    if np.any(uh < 0):
        n_neg = int(np.sum(uh < 0))
        logging.warning("UH has %d negative ordinates; clipping to zero.", n_neg)
        uh = np.clip(uh, 0.0, None)

    # --- Volume contract: a UH integrates to one unit (1 mm) of runoff over the area ---
    # 1 mm over area_km2 = area_km2 * 1e6 m2 * 0.001 m = area_km2 * 1000 m3
    expected_volume_m3 = area_km2 * 1000.0
    uh_volume_m3 = float(uh.sum() * dt_hours * 3600.0)
    ratio = uh_volume_m3 / expected_volume_m3 if expected_volume_m3 else np.nan
    logging.info("UH volume ratio (actual/expected for 1 mm): %.4f", ratio)
    if not np.isclose(ratio, 1.0, atol=0.02):
        logging.warning(
            "UH does not conserve unit volume (ratio %.3f). "
            "Rescaling ordinates so the UH represents exactly 1 mm of runoff.", ratio
        )
        uh = uh / ratio

    # --- Discrete convolution: Q_n = sum_m P_m * U_(n-m+1) ---
    direct_runoff = np.convolve(excess, uh, mode="full")
    logging.info(
        "Convolved %d excess pulses with %d UH ordinates -> %d DRH ordinates",
        excess.size, uh.size, direct_runoff.size,
    )

    # --- Volume balance check on the result ---
    runoff_volume_m3 = float(direct_runoff.sum() * dt_hours * 3600.0)
    input_volume_m3 = float(excess.sum() * expected_volume_m3 / 1.0)  # mm -> m3 over area
    balance = runoff_volume_m3 / input_volume_m3 if input_volume_m3 else np.nan
    logging.info(
        "Volume balance: runoff %.0f m3 vs input %.0f m3 (ratio %.4f)",
        runoff_volume_m3, input_volume_m3, balance,
    )

    total_flow = direct_runoff + baseflow_cms
    times = np.arange(direct_runoff.size) * dt_hours
    result = pd.DataFrame(
        {"direct_runoff_cms": direct_runoff, "total_flow_cms": total_flow},
        index=pd.Index(times, name="time_hours"),
    )
    logging.info("Peak total flow %.2f m3/s at t = %.2f h",
                 result["total_flow_cms"].max(), result["total_flow_cms"].idxmax())
    return result


if __name__ == "__main__":
    # Demo: a 3-pulse storm over a 25 km2 basin with a simple triangular UH
    dt = 1.0                     # hours
    basin_area = 25.0            # km2
    excess = np.array([5.0, 12.0, 6.0])   # mm per hour

    # Triangular UH (Tp = 3 h, Tb = 8 h) then normalised inside the function
    raw_uh = np.array([0, 2, 4, 6, 4.5, 3, 1.5, 0.5, 0.0])
    hydrograph = convolve_unit_hydrograph(
        excess_pulses_mm=excess,
        uh_ordinates=raw_uh,
        dt_hours=dt,
        area_km2=basin_area,
        baseflow_cms=1.5,
    )
    print(hydrograph.round(2))

The volume-ratio rescale is the safety net: it forces the UH to represent exactly one unit of runoff over the stated area regardless of how the raw ordinates were tabulated, so a mis-scaled synthetic UH cannot silently inflate the flood peak. np.convolve in full mode implements the discrete equation directly; scipy.signal.convolve(..., method="fft") is faster for very long records but returns the same ordinates.

Validation

Confirm four properties before trusting a convolved hydrograph:

Volume conservation. The volume under the direct-runoff hydrograph must equal the total excess depth times the basin area. In consistent units, sum(Q) * dt * 3600 (m³) should equal sum(P_mm) / 1000 * area_km2 * 1e6. A mismatch larger than one or two percent means the UH volume was wrong; the production function logs this ratio explicitly.

Ordinate count. The DRH must have exactly M + J − 1 ordinates. A shorter array means the convolution was truncated (mode="same" instead of mode="full"); a longer one means duplicate padding.

Peak timing. The DRH peak should fall after the most intense excess pulse by roughly the basin lag. A peak that coincides with the first pulse indicates the UH lost its rising limb.

Non-negativity. Every ordinate must be non-negative. Negative values point to a defective UH, usually from an unstable duration conversion — see the S-curve discussion above.

For calibration against observed storms, compare the simulated hydrograph to a gauge record using formal goodness-of-fit metrics rather than eyeballing the peak; the model calibration & objective functions topic covers Nash-Sutcliffe efficiency, KGE, and the volume/peak error statistics that quantify a UH fit.

Common Failure Modes

Duration mismatch. The most frequent and most damaging error. The excess pulse duration must equal the UH duration. Convolving a 1-hour UH with 15-minute excess inflates volume by a factor of four because each 15-minute pulse triggers a full 1-hour-unit response. Always S-curve the UH to the hyetograph time step first, and log both durations before convolving.

Non-conservation of volume. If the UH ordinates were built with the wrong cell area, time step, or unit-depth basis, the area under the UH is not one unit and every runoff pulse is scaled by the same wrong factor. Normalize the UH so sum(U) * dt equals one unit of depth over the area, and verify with a volume-balance check on the output.

Negative ordinates. These are physically impossible and usually arise from duration conversion. Subtracting two S-curve values that ring around their plateau produces negative differences; a UH derived by deconvolving a noisy observed storm can also go negative. Clip to zero and re-check the volume, or smooth the S-curve before differencing.

Ragged S-curve plateau. If the S-curve oscillates instead of settling flat at the supply rate, the UH ordinates are internally inconsistent with the stated duration. Any duration conversion built on that S-curve will propagate the ringing into the converted UH. Fix the source UH before converting.

Non-uniform excess assumption. The UH assumes excess rainfall is spatially uniform. On large or elongated basins a single storm is rarely uniform, and the lumped UH will misplace the peak. Subdivide the basin and route each sub-basin’s hydrograph, or move to a distributed model.

When to Use Unit Hydrographs vs Alternatives

Unit hydrograph convolution is the right tool for event-based design flood estimation on basins small enough to treat as spatially lumped, where you need a defensible peak discharge and hydrograph shape from limited data. It is fast, transparent, and standard in regulatory hydrology.

Reach for alternatives when its assumptions break:

• For the loss model that produces the excess input, and for composite runoff from mixed land cover, use curve number runoff estimation upstream of any UH.
• When you need a single-parameter synthetic UH for an ungauged basin, deriving SCS dimensionless unit hydrographs in Python is the most economical starting point.
• When choosing between the two dominant synthetic transforms — regional-coefficient Snyder versus GIS-distributed Clark — the trade-offs are laid out in Snyder vs Clark unit hydrograph: when to use each.
• For continuous simulation, spatially variable rainfall, or strongly nonlinear basins, move to a full hydrologic model such as those configured through HEC-HMS Python automation or the SWMM model integration with PySWMM, which handle routing, storage, and multi-event water balance that a single convolution cannot.

Related Topics

• Deriving SCS Dimensionless Unit Hydrographs in Python — the single-parameter NRCS method, peak-factor variants, and the t/Tp table implemented in numpy
• Snyder vs Clark Unit Hydrograph: When to Use Each — a decision guide contrasting regional-coefficient and time-area synthetic methods
• Curve Number Runoff Estimation — the loss model that produces the excess-rainfall input to convolution
• Model Calibration & Objective Functions — Nash-Sutcliffe, KGE, and error statistics for validating a UH against gauge data
• Rainfall-Runoff Modeling & Hydrologic Simulation — the parent overview linking losses, transforms, routing, and calibration