Watershed Modeling Python for hydrology automation

Green-Ampt vs Curve Number: Choosing a Loss Method

The loss method decides how much of a storm becomes runoff, and the two in common use answer that question from different information. The curve number method reads cumulative depth; Green-Ampt reads instantaneous intensity against a decaying infiltration capacity. Where those two views coincide, the methods agree. Where they diverge — short, intense storms on permeable ground — they can differ by a factor of two. This guide is part of the SCS curve number runoff estimation topic within rainfall-runoff modeling and hydrologic simulation.

Prerequisites

  • A hyetograph at the model time step — see precipitation forcing and design storms.
  • For the curve number: a composite CN for the subbasin.
  • For Green-Ampt: saturated conductivity, suction head, porosity and initial moisture content.
  • numpy, pandas.

Core Technique: Two Different Controlling Variables

What Each Method Is Watching A hyetograph with a short intense burst early and a long low-intensity tail. The curve number method produces no runoff until cumulative depth passes the initial abstraction, missing the early burst. Green-Ampt produces runoff during the burst, because intensity exceeds infiltration capacity there. rainfall intensity a short intense burst, then a long low-intensity tail runoff generated Green-Ampt: runoff during the burst curve number: nothing until Ia is met Same storm, same total depth, runoff generated at different times and in different amounts.

Annotated Code Example

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 curve_number_losses(hyeto_mm: pd.Series, cn: float,
                        ia_ratio: float = 0.2) -> pd.DataFrame:
    """
    SCS curve number losses, applied incrementally.

    The equation is defined on CUMULATIVE depth, so the incremental excess is
    the difference of the cumulative result — not the equation applied to each
    increment, which would be wrong and is a common implementation error.
    """
    S = 25400.0 / cn - 254.0            # potential retention, mm
    Ia = ia_ratio * S
    P = hyeto_mm.cumsum()

    Q = np.where(P > Ia, (P - Ia) ** 2 / (P - Ia + S), 0.0)
    excess = np.diff(np.concatenate([[0.0], Q]))

    log.info("CN %.0f → S = %.1f mm, Ia = %.1f mm", cn, S, Ia)
    log.info("Total rainfall %.1f mm → excess %.1f mm (runoff coefficient %.3f)",
             float(P.iloc[-1]), float(Q[-1]), float(Q[-1] / max(P.iloc[-1], 1e-9)))
    return pd.DataFrame({"rain_mm": hyeto_mm.to_numpy(),
                         "excess_mm": excess}, index=hyeto_mm.index)


def green_ampt_losses(
    hyeto_mm: pd.Series,
    timestep_min: float,
    ks_mm_hr: float,
    suction_mm: float,
    porosity: float,
    initial_moisture: float,
) -> pd.DataFrame:
    """
    Green-Ampt infiltration, stepped through the storm.

    Parameters
    ----------
    ks_mm_hr         : Saturated hydraulic conductivity.
    suction_mm       : Wetting-front suction head.
    porosity         : Soil porosity (volume fraction).
    initial_moisture : Initial volumetric moisture content.
    """
    dtheta = max(porosity - initial_moisture, 1e-4)
    dt_hr = timestep_min / 60.0

    F = 0.0                     # cumulative infiltration, mm
    excess, infil_rate = [], []

    for rain in hyeto_mm.to_numpy(dtype=float):
        i = rain / dt_hr if dt_hr > 0 else 0.0     # intensity, mm/hr

        # Infiltration CAPACITY at the current wetting-front depth. Early in
        # a storm F is small, so capacity is very large and all rain infiltrates.
        f_cap = ks_mm_hr * (1.0 + suction_mm * dtheta / F) if F > 0 else np.inf

        if i <= f_cap:
            # Supply-limited: everything infiltrates, no runoff whatever the
            # cumulative depth. This is where the two methods diverge most.
            f = i
        else:
            # Capacity-limited: ponding, and the excess becomes runoff.
            f = f_cap

        infiltrated = min(rain, f * dt_hr)
        F += infiltrated
        excess.append(max(0.0, rain - infiltrated))
        infil_rate.append(f if np.isfinite(f) else i)

    df = pd.DataFrame({"rain_mm": hyeto_mm.to_numpy(),
                       "excess_mm": excess,
                       "infil_capacity_mm_hr": infil_rate},
                      index=hyeto_mm.index)
    total_rain, total_excess = float(df.rain_mm.sum()), float(df.excess_mm.sum())
    log.info("Green-Ampt: Ks %.1f mm/hr, ψ %.0f mm, Δθ %.3f",
             ks_mm_hr, suction_mm, dtheta)
    log.info("Total rainfall %.1f mm → excess %.1f mm (runoff coefficient %.3f)",
             total_rain, total_excess, total_excess / max(total_rain, 1e-9))
    first = df.index[df.excess_mm > 0]
    if len(first):
        log.info("Runoff begins at minute %s, after %.1f mm of rain",
                 first[0], float(df.rain_mm.loc[:first[0]].sum()))
    return df


def compare_methods(hyeto_mm: pd.Series, timestep_min: float, cn: float,
                    **ga_params) -> pd.DataFrame:
    """Run both and report where and how much they differ."""
    cnr = curve_number_losses(hyeto_mm, cn)
    gar = green_ampt_losses(hyeto_mm, timestep_min, **ga_params)

    out = pd.DataFrame({
        "rain_mm": hyeto_mm,
        "excess_cn_mm": cnr["excess_mm"],
        "excess_ga_mm": gar["excess_mm"],
    })
    tot_cn, tot_ga = out.excess_cn_mm.sum(), out.excess_ga_mm.sum()
    log.info("Total excess — CN %.1f mm, Green-Ampt %.1f mm (%+.1f %%)",
             tot_cn, tot_ga, 100.0 * (tot_ga - tot_cn) / max(tot_cn, 1e-9))

    peak_cn = out.excess_cn_mm.max() / (timestep_min / 60.0)
    peak_ga = out.excess_ga_mm.max() / (timestep_min / 60.0)
    log.info("Peak excess intensity — CN %.1f mm/hr, Green-Ampt %.1f mm/hr",
             peak_cn, peak_ga)
    if abs(tot_ga - tot_cn) / max(tot_cn, 1e-9) > 0.25:
        log.warning("The two methods differ by more than 25 %% on this storm — "
                    "this is the intensity-sensitive regime, so the choice of "
                    "method is a modelling decision, not a detail")
    return out


# --- Example usage ---
# out = compare_methods(hyeto, timestep_min=15, cn=78,
#                       ks_mm_hr=10.9, suction_mm=110.0,
#                       porosity=0.463, initial_moisture=0.30)

Parameter Reference

Curve number Green-Ampt
Parameters 1 (CN), plus an Ia ratio 4 (Ks, ψ, porosity, θᵢ)
Controlling variable Cumulative depth Instantaneous intensity
Responds to storm shape No Yes
Antecedent moisture Via an AMC class shift Via θᵢ, continuously
Source of parameters Land cover × soil group lookup Soil texture class tables
Parameter uncertainty ±5 CN units is typical Ks spans an order of magnitude within a class
Continuous simulation Needs external moisture accounting Natural, with recovery between storms
Regulatory acceptance Very wide Wide, but less universal

Worked Example: Where They Split

The same 152 mm design storm on the same subbasin, run at three storm durations:

Duration Peak intensity CN excess Green-Ampt excess Difference
24 h 45 mm/h 78.4 mm 74.1 mm −5 %
6 h 96 mm/h 78.4 mm 89.2 mm +14 %
1 h 220 mm/h 78.4 mm 118.6 mm +51 %

The curve number column is constant by construction: same depth, same answer. Green-Ampt’s rises steeply as the same water arrives faster, because infiltration capacity is exceeded for more of the storm. On a permeable soil under a short intense storm, the two methods are answering different questions and the difference is not a calibration problem.

The Same Depth, Delivered Faster With a fixed 152 millimetre total, the curve number excess stays at 78.4 millimetres for every duration. Green-Ampt excess rises from 74.1 at 24 hours to 89.2 at 6 hours and 118.6 at 1 hour, because higher intensity exceeds infiltration capacity for more of the storm. 24 h 78.4 vs 74.1 6 h 78.4 vs 89.2 1 h 78.4 vs 118.6 curve number Green-Ampt excess mm 152 mm total in every case

Green-Ampt’s four parameters are not equally uncertain. Sweeping each across its published range for one soil class shows that conductivity dominates and the other three barely register.

One Parameter Carries Almost All the Uncertainty Sweeping saturated conductivity across its published range for a silt loam moves the computed runoff by 62 percent. Initial moisture moves it by 19 percent, porosity by 6 and suction head by 4. saturated conductivity 62 % initial moisture 19 % porosity 6 % suction head 4 % change in computed runoff across each parameter's published range Effort spent narrowing the last three is effort not spent on the first.

Gotchas and Edge Cases

  • Applying the CN equation per increment. The equation is defined on cumulative depth; differencing the cumulative result is the correct incremental excess, and applying it per increment overstates runoff badly.
  • Green-Ampt with no ponding logic. Rain that arrives below capacity must all infiltrate, whatever the cumulative total. Implementations that compare cumulative depth against capacity produce curve-number behaviour under a Green-Ampt name.
  • Ks taken from the class midpoint. Published conductivity ranges span an order of magnitude within a texture class; the midpoint is a guess with a very wide interval, and the runoff estimate inherits it.
  • Initial moisture left at a default. It is the antecedent-condition parameter, and it moves the answer as much as an AMC class shift does in the curve number method.
  • Comparing on a 24-hour design storm only. The two methods agree there and diverge on short storms. Comparing on the case where they agree tells you nothing about the case where you need them.
  • Switching methods mid-calibration. The parameters are not interchangeable, and a model calibrated with one loss method is not calibrated with the other.