Watershed Modeling Python for hydrology automation

Precipitation Forcing and Design Storms

Every rainfall-runoff model is a function whose most important argument is rainfall, and more calibration effort is wasted compensating for bad forcing than on any other cause. As part of the rainfall-runoff modeling and hydrologic simulation workflow, this stage decides what the model is actually being asked to reproduce — and a 15 % error in basin-average depth cannot be recovered by adjusting a curve number or a unit hydrograph, only disguised by them.

There are two distinct jobs here that are often conflated. Observed forcing reconstructs what fell during a real event, for calibration and validation. Design forcing synthesises a storm that has never happened, for a design flood of a stated recurrence interval. They use different data, different assumptions and different validation, and a pipeline that treats them as one thing produces a model calibrated on reality and applied to a fiction that does not match it.

Prerequisites and Environment Setup

bash 
conda create -n forcing python=3.11
conda activate forcing
conda install -c conda-forge xarray=2024.2 rioxarray=0.15 pandas=2.2 \
    geopandas=0.14 requests netcdf4 numpy scipy matplotlib
Input Requirement Notes
Catchment boundary Polygon in a projected CRS Areal averaging and reduction factors both need real area
Time of concentration Estimated for the catchment Sets the minimum storm duration
Depth-duration-frequency source NOAA Atlas 14, or the regional equivalent Point depths by return period and duration
Model time step Known before the hyetograph is built A hyetograph must be resampled to it, not interpolated afterwards
Gauge or gridded record For calibration events only Hourly or finer; daily totals cannot drive an event model

The time step deserves a note. A hyetograph built at 15-minute resolution and then handed to a model running at 1-hour steps loses its peak intensity to averaging, and the modelled peak falls with it. Build the hyetograph at the step the model will run.

Mechanics: Where Rainfall Data Comes From

Four families of precipitation data reach hydrologic models, and they differ in what they actually measure.

What Each Precipitation Product Actually Measures Gauges measure depth at a point with high accuracy but no spatial coverage. Radar measures reflectivity aloft, converted to rate, with high resolution and systematic bias. Gauge-adjusted radar combines both. Interpolated grids such as PRISM measure long-term climatology well but smooth individual storms. product measures resolution characteristic bias rain gauge depth at one point 1 min, zero area under-catch in wind raw radar reflectivity aloft 5 min, 1 km Z-R error, tens of % gauge-adjusted radar pattern + depth 1 h, 4 km residual, single digits interpolated grid smoothed depth field 1 day, 800 m peaks smoothed away A daily interpolated grid cannot drive an event model however fine its spatial resolution looks.

The distinction that matters most is between products that measure depth and products that infer it. A gauge measures depth. Radar measures reflectivity in a volume some hundreds of metres above the ground and converts it to a rate through a power law whose coefficients vary by storm type; the resulting bias is systematic, not random, and it does not average out over an event.

Point depths, basin depths, and the factor between them

A depth-duration-frequency product gives the depth expected at a point. A storm covering 400 km² does not deliver its point-maximum depth everywhere, so basin-average depth is lower — and the gap grows with basin area and shrinks with duration, because long storms are more spatially uniform than short convective ones.

Areal Reduction Grows with Area and Shrinks with Duration Four curves of areal reduction factor against basin area from 10 to 1000 square kilometres. A fifteen-minute storm falls to about 0.55 at 1000 square kilometres, a one-hour storm to 0.7, a six-hour storm to 0.85 and a twenty-four-hour storm to about 0.93. 10 50 150 400 1 000 basin area (km², log scale) 0.5 0.6 0.8 0.9 1.0 reduction factor 15 min 1 h 6 h 24 h below ~25 km² all four sit above 0.95 Skipping the reduction on a 400 km² basin overstates a one-hour design depth by roughly a quarter.

Step-by-Step Workflow

  1. Fix the return period and duration. The return period comes from the design standard. The duration must be at least the catchment’s time of concentration — a storm shorter than the travel time never engages the whole basin, and the resulting peak is not the design peak. See estimating time of concentration with the NRCS velocity method.
  2. Retrieve the depth-duration-frequency depths at the catchment centroid for every sub-duration you will need — typically 5 min through the full storm duration.
  3. Apply the areal reduction factor appropriate to the basin area and the storm duration.
  4. Distribute the depth in time. Either an SCS type curve, which encodes a regional storm shape, or the alternating-block method, which constructs a hyetograph consistent with the whole depth-duration-frequency curve.
  5. Resample to the model time step by aggregation, not interpolation.
  6. Verify. The hyetograph must sum to the design depth, and each sub-duration window must still match its own depth.

The two distribution methods

SCS type curves are fixed dimensionless mass curves — Type I, IA, II and III — each representing a regional storm shape over 24 hours. They are simple, defensible where the regional type is agreed, and blunt: every 24-hour storm in a region gets the same shape regardless of return period.

The alternating-block method builds a hyetograph directly from the depth-duration-frequency curve. Compute the depth at each successive duration, difference them to get incremental depths, then place the largest increment centrally with subsequent increments alternating either side. The result satisfies every sub-duration depth simultaneously, which is why it is preferred for durations where a type curve does not exist.

How the Alternating-Block Hyetograph Is Assembled Twelve 15-minute increments derived by differencing cumulative depth-duration-frequency depths. The largest increment is placed at position 6, the second largest at position 7, the third at position 5, and so on alternating outward, producing a single-peaked hyetograph. 1st 2nd 3rd 4th 5th 6th 0 1.5 h 3 h time — twelve 15-minute increments depth per step Placement order is what makes every centred window match its own design depth simultaneously.

Production-Ready Code

The function below builds an alternating-block hyetograph from a depth-duration-frequency table, applies an areal reduction factor and returns a pandas series at the model time step.

python 
import logging

import numpy as np
import pandas as pd

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


def areal_reduction_factor(area_km2: float, duration_min: float) -> float:
    """
    A smooth areal reduction curve of the TP-29 family.

    Returns 1.0 for very small basins and decays with area, more slowly for
    long durations. Regional guidance should override this where it exists —
    this is a defensible default, not a substitute for a published curve.
    """
    if area_km2 <= 2.5:
        return 1.0
    # Longer storms are more spatially uniform, so they decay more slowly.
    decay = 0.048 * np.exp(-0.0032 * duration_min) + 0.006
    arf = float(np.exp(-decay * np.power(area_km2, 0.55)))
    return max(0.4, min(1.0, arf))


def alternating_block_hyetograph(
    ddf_depths_mm: dict[int, float],
    duration_min: int,
    timestep_min: int,
    area_km2: float | None = None,
) -> pd.Series:
    """
    Build a design hyetograph consistent with a whole depth-duration-frequency curve.

    Parameters
    ----------
    ddf_depths_mm : Cumulative point depth by duration in minutes, e.g.
                    {5: 12.1, 10: 18.4, 15: 23.0, 30: 32.6, 60: 41.9, 120: 50.3}.
    duration_min  : Total storm duration; must be a key in ddf_depths_mm.
    timestep_min  : Model time step; must divide duration_min exactly.
    area_km2      : Basin area. When given, an areal reduction factor is applied.

    Returns
    -------
    A pandas Series of incremental depth (mm) indexed by minutes from storm start.
    """
    if duration_min % timestep_min:
        raise ValueError("duration must be a whole number of time steps")
    if duration_min not in ddf_depths_mm:
        raise ValueError(f"no depth-duration-frequency depth for {duration_min} min")

    n_steps = duration_min // timestep_min
    durations = [timestep_min * (i + 1) for i in range(n_steps)]

    # --- Interpolate the cumulative depth curve onto our step boundaries.
    # Interpolating in log-log space keeps the curve monotonic and smooth,
    # which matters because differencing it must not produce negatives. ---
    known_d = np.array(sorted(ddf_depths_mm))
    known_p = np.array([ddf_depths_mm[d] for d in known_d])
    cumulative = np.exp(np.interp(np.log(durations), np.log(known_d), np.log(known_p)))

    arf = 1.0
    if area_km2 is not None:
        arf = areal_reduction_factor(area_km2, duration_min)
        cumulative = cumulative * arf
        log.info("Areal reduction factor for %.1f km² over %d min: %.3f",
                 area_km2, duration_min, arf)

    # --- Incremental depths: the extra rain each successive window adds ---
    increments = np.diff(np.concatenate([[0.0], cumulative]))
    if (increments < 0).any():
        log.warning("Non-monotonic depth-duration curve — %d negative increment(s) "
                    "clipped to zero; check the source table", int((increments < 0).sum()))
        increments = np.clip(increments, 0.0, None)

    # --- Alternating placement: largest centrally, then outward either side ---
    order = np.argsort(increments)[::-1]        # largest first
    slots = [0] * n_steps
    centre = n_steps // 2
    left, right = centre - 1, centre
    for rank, idx in enumerate(order):
        if rank == 0:
            slots[centre] = increments[idx]
            continue
        if rank % 2 == 1:
            slots[right + 1] = increments[idx]
            right += 1
        else:
            slots[left] = increments[idx]
            left -= 1

    index = [timestep_min * i for i in range(n_steps)]
    series = pd.Series(slots, index=index, name="depth_mm")
    series.index.name = "minutes_from_start"

    total = float(series.sum())
    log.info("Hyetograph: %d steps of %d min, total %.1f mm (design %.1f mm, ARF %.3f)",
             n_steps, timestep_min, total, ddf_depths_mm[duration_min], arf)
    log.info("Peak intensity %.1f mm/h at minute %d",
             series.max() * 60 / timestep_min, int(series.idxmax()))
    return series


# --- Example usage: a 100-year, 3-hour storm on a 180 km² basin ---
# ddf = {5: 14.2, 10: 21.6, 15: 27.1, 30: 38.4, 60: 49.8, 120: 60.2, 180: 66.5}
# hyeto = alternating_block_hyetograph(
#     ddf_depths_mm=ddf, duration_min=180, timestep_min=15, area_km2=180.0
# )
# print(hyeto.round(2))

Validation Protocol

Forcing errors are silent — the model runs, the hydrograph looks like a hydrograph — so the checks have to be explicit.

  • Total depth. The hyetograph must sum to the design depth after areal reduction. An error here is arithmetic and inexcusable.
  • Every sub-duration window. For an alternating-block hyetograph, the maximum rolling sum over each sub-duration must equal that duration’s depth. This is the check that proves the construction, and it fails immediately if the increments were placed in the wrong order.
  • Peak intensity against the source. The maximum intensity should equal the shortest-duration depth divided by the time step. If it is lower, the hyetograph was interpolated to the model step rather than aggregated.
  • Basin-average depth against gauges, for observed events. Where the basin holds gauges, compare the basin-average depth from the gridded product against the gauge mean. A systematic offset is a product bias and must be corrected before calibration, not during it.
  • Volume closure against runoff. For a calibration event, the runoff volume cannot exceed the rainfall volume. A model that needs a runoff coefficient above 1.0 to match the observed hydrograph has a forcing error, not a parameter problem.

Common Failure Modes and Optimization

  • Point depths applied to a large basin. The most common single error, and it always inflates the design flood. On a 400 km² basin with a one-hour storm the overstatement is around 25 %.
  • Storm shorter than the time of concentration. The whole basin never contributes simultaneously, so the computed peak is not the design peak. Always check that duration exceeds the travel time.
  • Hyetograph interpolated to the model step. Interpolation smooths the peak; aggregation preserves the volume and the correct intensity. Build at the model step from the start.
  • Daily products driving event models. A daily total spread uniformly over 24 hours has a peak intensity an order of magnitude below reality, and no loss method compensates for that.
  • Raw radar for calibration. Reflectivity converted through a fixed power law carries a storm-dependent bias. Use a gauge-adjusted product, or adjust it yourself against the gauges in the basin.
  • Time zone and interval convention. A gridded product stamped at the end of its accumulation interval, read as though stamped at the start, shifts the whole hyetograph by one step. This shows up as a systematic timing error that calibration then “fixes” by distorting the transform.

When to Use This vs. Alternatives

Use a design storm when the question is about a recurrence interval — a culvert size, a floodplain extent, a detention volume. There is no observed event with a stated return period, so a synthetic one is the only option.

Use observed forcing when the question is about model skill. Calibration and validation both require events that actually happened, with observed flow to compare against. See model calibration objective functions for what to do with them.

Use continuous forcing when the question involves antecedent conditions, yield or low flow. A single design storm carries no memory, and the assumed antecedent moisture is doing far more work than most reports admit.

Skip the rainfall model entirely when a gauged flood-frequency analysis will answer the question directly. If the site has a long gauge record, fitting a frequency curve to the observed peaks is more defensible than a rainfall-runoff chain with a design storm at one end — see flood frequency and streamflow statistics.

Frequently Asked Questions

What is an areal reduction factor and when must I apply one?

Depth-duration-frequency values from products like NOAA Atlas 14 are point depths — the depth expected at a single location. A storm never delivers its point-maximum depth uniformly across a large basin, so applying a point depth to a whole catchment overstates the volume. The areal reduction factor scales the point depth down as a function of basin area and storm duration. Below about 25 km² it is close to 1.0 and often ignored; at 500 km² for a one-hour storm it can fall below 0.7, which is far too large to neglect.

Why does an alternating-block hyetograph put the peak in the middle?

The method is constructed so that every sub-duration window centred on the peak matches its depth-duration-frequency depth. Placing the largest increment centrally, then alternating the next largest to either side, is what makes the 15-minute, 1-hour and 6-hour windows all simultaneously correct. It is a mathematical consequence of the construction, not a claim that real storms peak in the middle.

Should I use radar or gauge precipitation to calibrate a model?

Use gauge data where the gauge network resolves the storm and the basin is small, because gauges measure depth directly. Use gauge-adjusted radar for larger basins and convective storms, where a single gauge cannot represent the spatial pattern. Raw radar reflectivity converted through a fixed power law should not be used for calibration at all — its bias is often tens of percent and varies by storm type.