Rainfall-Runoff Modeling & Hydrologic Simulation in Python

Rainfall-runoff modeling is the engineering discipline of turning a precipitation record into a simulated streamflow hydrograph through a reproducible computational pipeline. The problem is deceptively simple to state and hard to solve well: given how much rain fell, where, and how fast, predict the discharge that will pass a downstream point as a function of time. This section treats that problem as a Python engineering workflow rather than a black-box modeling exercise, chaining basin parameterization, loss accounting, hydrograph transforms, baseflow separation, and channel routing into scripts that run identically on a laptop and in the cloud. The five topics that make up this section each own one link in that chain: HEC-HMS Python automation for driving the industry-standard engine, SWMM model integration with pyswmm for urban and stormwater systems, unit hydrograph methods for the rainfall-to-runoff transform, curve number runoff estimation for the loss calculation, and model calibration with objective functions for tuning the whole system against observations. Every model here consumes upstream products: the subbasin geometry and routing topology produced by watershed delineation and catchment synchronization, and the conditioned terrain and derived rasters from hydrology data preparation and DEM processing. A hydrologic model is only as trustworthy as the basin parameters it inherits from those two stages.

Foundational Concepts

The Hydrologic Cycle as a Computational Abstraction

Every rainfall-runoff model is a discretized bookkeeping of the same physical sequence. Precipitation arrives at the land surface, where a fraction is intercepted by canopy and depression storage and a further fraction infiltrates into the soil. Whatever remains after these losses is excess rainfall — the water that actually becomes surface runoff. That excess is routed off the hillslopes and concentrated into a channel, a process the model represents as a transform from a depth-over-time signal into a discharge-over-time signal (the direct runoff hydrograph). Slower subsurface pathways contribute baseflow, a smooth recession that is superimposed on the direct runoff. Finally, the combined flow travels down the channel network, where routing attenuates and delays the wave as it moves from reach to reach.

Expressed as a pipeline, the chain is: precipitation, then interception and infiltration losses, then excess rainfall, then transform to a hydrograph, then baseflow addition, then channel routing, then the simulated hydrograph at the outlet. Each arrow in that chain is a numerical operator with parameters, and the entire engineering task reduces to estimating those parameters defensibly and evaluating the operators in the right order on a consistent timestep.

Lumped, Semi-Distributed, and Distributed Models

Models differ chiefly in how finely they resolve space:

• Lumped models treat an entire watershed as a single unit with one set of parameters (one curve number, one time of concentration). They are fast, easy to calibrate, and appropriate when the basin is small and reasonably homogeneous. A single unit hydrograph applied to a basin-average excess rainfall is the canonical lumped approach.
• Semi-distributed models partition the watershed into subbasins, each lumped internally but connected through a channel routing network. This is the workhorse structure of HEC-HMS: parameters vary between subbasins, and the routing reaches carry each subbasin’s hydrograph to a common outlet. The subbasin partition comes directly from basin partitioning strategies.
• Distributed models solve the governing equations on a grid, letting rainfall, soils, and land cover vary cell by cell. They demand the most data and computation and are justified when spatial heterogeneity is the point — for example, gridded surface routing in a fully distributed rainfall-runoff engine.

Choosing the right spatial resolution is a data-availability decision as much as an accuracy one. A distributed model fed with poorly characterized soils will not outperform a well-calibrated semi-distributed one.

Event versus Continuous Simulation

An event model simulates a single storm. It assumes an antecedent moisture condition, applies a loss method for the duration of the event, and stops. It is the right tool for design-storm peak flows and flood hydrographs, and it is cheap to calibrate because it ignores what happens between storms. A continuous model runs across months or years, carrying soil-moisture and groundwater storage forward from one storm to the next. Continuous simulation is required when antecedent conditions, seasonal baseflow, or long-term water balance drive the answer, but it needs continuous forcing (precipitation and evapotranspiration) and many more parameters. The distinction shapes every downstream decision, from which loss method is admissible to how the model is validated.

Model and Method Landscape

The methods in this section are not competitors so much as interchangeable operators that occupy specific slots in the modeling chain. A loss method, a transform, and a routing method combine to form a complete model; HEC-HMS and SWMM are engines that bundle those operators together with a solver and a user interface.

The two loss methods illustrate the trade-off cleanly. The curve number method is empirical, needs only a composite CN derived from land cover and hydrologic soil group, and is the regulatory default for event modeling — but it has no explicit time dimension and behaves poorly for very small storms. Green-Ampt is physically based and resolves infiltration through the event, at the cost of soil-hydraulic parameters that are hard to measure. Similarly, the choice between a Snyder and a Clark unit hydrograph, or the derivation of an SCS dimensionless unit hydrograph, determines the shape of the transform without changing the loss accounting that precedes it. Kinematic-wave routing is preferred where channel slopes are steep and backwater effects are negligible; Muskingum and Muskingum-Cunge remain the pragmatic default for gentler reaches because they need only two parameters that calibrate readily.

Production Python Architecture

A production rainfall-runoff pipeline has three stages: derive basin parameters from geospatial data, drive a model engine, and post-process the resulting hydrograph. The first stage uses rasterio and geopandas to reduce a subbasin polygon and the rasters underneath it to scalar parameters — drainage area, composite curve number, time of concentration, and lag. The NRCS velocity method for time of concentration and composite curve numbers from land-use rasters each get a dedicated treatment; the block below shows the transform, loss, and hydrograph assembly as a self-contained, runnable model that depends only on numpy and pandas.

The engine stage can either call this in-process transform or drive an external solver. For urban systems you would hand the subcatchment and conduit network to pyswmm and step the simulation, as covered in running and post-processing SWMM simulations with pyswmm. For HEC-HMS you write a control script and invoke the engine through its Jython scripting API from a Python subprocess, then read the resulting DSS time series back into pandas.

import logging
import numpy as np
import pandas as pd

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s [%(levelname)s] %(name)s — %(message)s",
)
log = logging.getLogger("rainfall_runoff")


def scs_excess_rainfall(rainfall_mm: np.ndarray, curve_number: float) -> np.ndarray:
    """
    Convert incremental rainfall depths to incremental excess (runoff) depths
    using the SCS curve number method.

    Parameters
    ----------
    rainfall_mm : np.ndarray
        Incremental rainfall depth per timestep, millimetres.
    curve_number : float
        Composite curve number (30–98) for the subbasin.

    Returns
    -------
    np.ndarray
        Incremental excess rainfall depth per timestep, millimetres.
    """
    if not 30.0 <= curve_number <= 98.0:
        raise ValueError(f"Curve number {curve_number} outside valid 30–98 range.")

    # Potential maximum retention S and initial abstraction Ia (Ia = 0.2 S)
    s = 25400.0 / curve_number - 254.0  # mm
    ia = 0.2 * s
    log.info("CN=%.1f -> S=%.1f mm, Ia=%.1f mm", curve_number, s, ia)

    cumulative_p = np.cumsum(rainfall_mm)
    # Cumulative runoff Q is zero until cumulative P exceeds Ia
    runoff_cum = np.where(
        cumulative_p > ia,
        (cumulative_p - ia) ** 2 / (cumulative_p - ia + s),
        0.0,
    )
    excess = np.diff(runoff_cum, prepend=0.0)
    log.info("Total rainfall %.1f mm -> total excess %.1f mm",
             cumulative_p[-1], runoff_cum[-1])
    return excess


def scs_unit_hydrograph(area_km2: float, tc_hr: float, dt_hr: float) -> np.ndarray:
    """
    Build an SCS dimensionless unit hydrograph scaled to a subbasin.

    The peak of the UH occurs at time-to-peak Tp = 0.6 * Tc + dt/2, and the
    peak discharge follows the standard SCS 484 relationship.
    """
    lag_hr = 0.6 * tc_hr
    tp_hr = dt_hr / 2.0 + lag_hr
    qp = 2.08 * area_km2 / tp_hr  # m3/s per mm of excess (SI 484 form)
    log.info("UH: Tp=%.2f h, peak=%.3f m3/s per mm excess", tp_hr, qp)

    # SCS dimensionless UH ordinates (t/Tp -> q/qp), interpolated to dt grid
    t_ratio = np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6,
                        1.8, 2.0, 2.4, 2.8, 3.2, 4.0, 5.0])
    q_ratio = np.array([0.0, 0.10, 0.31, 0.66, 0.93, 1.0, 0.93, 0.78, 0.56,
                        0.39, 0.28, 0.147, 0.077, 0.040, 0.011, 0.0])
    t_end = t_ratio[-1] * tp_hr
    t_grid = np.arange(0.0, t_end + dt_hr, dt_hr)
    uh = np.interp(t_grid / tp_hr, t_ratio, q_ratio) * qp
    return uh


def route_muskingum(inflow: np.ndarray, k_hr: float, x: float, dt_hr: float) -> np.ndarray:
    """
    Route an inflow hydrograph through a channel reach with the Muskingum method.

    Parameters
    ----------
    k_hr : float
        Travel time (storage) constant, hours.
    x : float
        Weighting factor (0.0–0.5); 0 gives maximum attenuation.
    """
    denom = k_hr - k_hr * x + 0.5 * dt_hr
    c0 = (0.5 * dt_hr - k_hr * x) / denom
    c1 = (0.5 * dt_hr + k_hr * x) / denom
    c2 = (k_hr - k_hr * x - 0.5 * dt_hr) / denom
    log.info("Muskingum coefficients C0=%.3f C1=%.3f C2=%.3f (sum=%.3f)",
             c0, c1, c2, c0 + c1 + c2)

    outflow = np.zeros_like(inflow)
    outflow[0] = inflow[0]
    for t in range(1, len(inflow)):
        outflow[t] = c0 * inflow[t] + c1 * inflow[t - 1] + c2 * outflow[t - 1]
    return np.clip(outflow, 0.0, None)


def simulate_hydrograph(
    rainfall_mm: np.ndarray,
    area_km2: float,
    curve_number: float,
    tc_hr: float,
    dt_hr: float,
    baseflow_m3s: float = 0.5,
    k_hr: float = 2.0,
    x: float = 0.2,
) -> pd.DataFrame:
    """
    End-to-end event simulation: losses, transform, baseflow, and routing.
    """
    log.info("Simulating basin: area=%.2f km2, CN=%.1f, Tc=%.2f h, dt=%.2f h",
             area_km2, curve_number, tc_hr, dt_hr)

    excess = scs_excess_rainfall(rainfall_mm, curve_number)
    uh = scs_unit_hydrograph(area_km2, tc_hr, dt_hr)

    # Direct runoff = convolution of excess rainfall with the unit hydrograph
    direct = np.convolve(excess, uh)[: len(excess)]
    total_inflow = direct + baseflow_m3s
    routed = route_muskingum(total_inflow, k_hr, x, dt_hr)

    peak = float(routed.max())
    t_peak = float(np.argmax(routed) * dt_hr)
    log.info("Simulated peak %.2f m3/s at t=%.2f h", peak, t_peak)

    time_hr = np.arange(len(routed)) * dt_hr
    return pd.DataFrame(
        {
            "time_hr": time_hr,
            "excess_mm": np.pad(excess, (0, len(routed) - len(excess))),
            "direct_m3s": direct,
            "routed_m3s": routed,
        }
    )


if __name__ == "__main__":
    # Synthetic 6-hour design storm on a 15-minute timestep
    dt = 0.25
    storm = np.array([2, 4, 8, 15, 22, 18, 12, 8, 5, 3, 2, 1] + [0] * 60,
                     dtype=float)
    result = simulate_hydrograph(
        rainfall_mm=storm,
        area_km2=42.5,
        curve_number=78.0,
        tc_hr=3.4,
        dt_hr=dt,
    )
    log.info("Hydrograph rows: %d, runoff volume %.0f m3",
             len(result), result["routed_m3s"].sum() * dt * 3600.0)

The function chain is deliberately transparent: each operator is testable in isolation, the Muskingum coefficients are logged so a reviewer can confirm they sum to one, and the water balance can be checked by comparing total excess depth against routed volume. When an external engine replaces the in-process transform, only simulate_hydrograph changes — the parameter-derivation and post-processing stages stay identical, which is what makes the pattern portable across HEC-HMS and SWMM.

Scaling and Cloud-Native Workflows

A single event simulation runs in milliseconds, but production hydrology rarely needs just one. Three scaling patterns cover most workloads.

Batch event runs. Design-storm analysis requires simulating a matrix of return periods and durations (the 2-, 10-, 25-, and 100-year storms at 6-, 12-, and 24-hour durations). Structure this as a parameter grid and map the simulation function over it. Because each run is independent, the batch is embarrassingly parallel; a concurrent.futures.ProcessPoolExecutor saturates a workstation, and the same callable submits to a distributed pool of compute nodes without modification.

Monte-Carlo and ensemble forcing. Uncertainty analysis perturbs the forcing (rainfall depth, temporal distribution) and the parameters (curve number, Tc, Muskingum K) across thousands of realizations to produce a distribution of peak flows rather than a single deterministic value. Ensemble numerical-weather-prediction forcing feeds the same machinery for probabilistic flood forecasting. Store the realization matrix in a pandas DataFrame or a Zarr array and stream results back as they complete.

import logging
import numpy as np
from concurrent.futures import ProcessPoolExecutor

log = logging.getLogger("ensemble")


def monte_carlo_peaks(base_storm: np.ndarray, n_real: int = 2000, seed: int = 42):
    """Sample CN and Tc, simulate, and return the peak-flow distribution."""
    rng = np.random.default_rng(seed)
    cn_samples = rng.normal(78.0, 4.0, n_real).clip(50, 95)
    tc_samples = rng.normal(3.4, 0.5, n_real).clip(0.5, 12.0)
    log.info("Launching %d Monte-Carlo realizations", n_real)

    def _one(args):
        cn, tc = args
        df = simulate_hydrograph(base_storm, 42.5, cn, tc, 0.25)
        return float(df["routed_m3s"].max())

    with ProcessPoolExecutor() as pool:
        peaks = np.fromiter(pool.map(_one, zip(cn_samples, tc_samples)),
                            dtype=float)
    log.info("Peak flow: median %.1f, 95th pct %.1f m3/s",
             np.median(peaks), np.percentile(peaks, 95))
    return peaks

Dask-parallel calibration. Automated calibration evaluates the model hundreds to thousands of times inside an optimizer, and those evaluations dominate wall-clock time. A dask.distributed scheduler spreads the objective-function evaluations across a pool of workers, letting global search algorithms such as differential evolution or the SCE-UA method explore the parameter space in parallel. The automated calibration workflow with SciPy and spotpy covers the optimizer wiring. Where the pipeline stitches together delineation and modeling for many basins, orchestration with Prefect or Airflow keeps the DEM-to-hydrograph chain reproducible, as described in pipeline orchestration from DEM to watershed.

Validation and QA/QC

A simulated hydrograph is a hypothesis until it is scored against observations. The reference data are almost always continuous discharge records from stream gauges, most commonly the USGS National Water Information System at https://waterdata.usgs.gov/, which serves instantaneous and daily values that align directly with model timesteps. Three complementary objective functions quantify agreement, and no single one is sufficient on its own:

• Nash-Sutcliffe efficiency (NSE) compares the model’s squared error against the variance of the observations. NSE of 1.0 is perfect, 0.0 means the model is no better than the observed mean, and negative values mean it is worse. Because it squares errors, NSE is dominated by high flows and rewards matching peaks and hydrograph shape.
• Kling-Gupta efficiency (KGE) decomposes performance into correlation, bias ratio, and variability ratio. This decomposition is diagnostically valuable: a mediocre KGE tells you not just that the model is wrong but whether the timing, the volume, or the flashiness is at fault.
• Percent bias (PBIAS) measures systematic volume error as a signed percentage. Positive PBIAS indicates under-prediction of total runoff, negative indicates over-prediction. It is the cleanest check on water-balance closure.

The computation of NSE and KGE in Python walks through the exact formulas and edge cases. In practice, calibrate against KGE or a weighted blend, then report NSE, KGE, and PBIAS together so a reviewer can judge timing and volume independently. Beyond aggregate metrics, always inspect the hydrograph visually: overlay simulated and observed discharge on the same axes and confirm that the rising limb, peak timing, peak magnitude, and recession all track. A model can score a respectable NSE while getting the peak timing badly wrong if the record is dominated by low flows.

Common Failure Modes

Water-balance closure error. Root cause: the total simulated runoff volume does not equal precipitation minus losses, usually from an inconsistent unit conversion between rainfall depth (mm), basin area (km²), and discharge (m³/s), or from a truncated hydrograph tail. Mitigation: assert that integrated routed volume equals excess-rainfall depth times basin area to within a small tolerance, and extend the simulation window until the recession limb reaches baseflow before integrating.

Unrealistic time of concentration. Root cause: Tc estimated from an incorrect longest-flow-path length or slope produces a peak that arrives too early or too late and a hydrograph that is too flashy or too damped. Mitigation: derive the flow path from the delineated network rather than a straight-line distance, cross-check the NRCS velocity-method Tc against an independent lag equation, and confirm the model timestep is at most one-fifth of Tc so the rising limb is resolved.

Over-parameterization and equifinality. Root cause: fitting too many free parameters against a single event lets many different parameter sets reproduce the same hydrograph, so the calibration is non-unique and does not generalize. Mitigation: fix physically constrained parameters (area, CN bounds) from data, calibrate only the few that are genuinely uncertain, and validate on withheld events. Global sensitivity analysis before calibration reveals which parameters actually matter.

Timestep instability. Root cause: a routing timestep that is too large relative to the reach travel time violates numerical stability, producing negative flows or oscillations, particularly in kinematic-wave and Muskingum routing. Mitigation: keep the Muskingum weighting factor and timestep within the stability envelope (verify the routing coefficients stay non-negative), and satisfy the Courant condition for kinematic-wave routing by shrinking the timestep or subdividing long reaches.

Related Topics

• HEC-HMS Python Automation — scripting the industry-standard engine for batch and reproducible runs
• SWMM Model Integration with pyswmm — driving and post-processing urban stormwater simulations from Python
• Unit Hydrograph Methods — the transform from excess rainfall to a direct runoff hydrograph
• Curve Number Runoff Estimation — SCS loss accounting and time-of-concentration estimation
• Model Calibration & Objective Functions — NSE, KGE, and automated parameter optimization
• Watershed Delineation & Catchment Synchronization — the source of subbasin geometry and routing topology that hydrologic models consume
• Hydrology Data Preparation & DEM Processing — the conditioned terrain and derived rasters that feed basin parameterization