SCS Curve Number Runoff Estimation in Python

The NRCS curve number method converts a rainfall depth into a runoff depth using a single empirical parameter that folds soil infiltration capacity, land cover, and antecedent moisture into one number between roughly 30 and 100. Within the Rainfall-Runoff Modeling & Hydrologic Simulation domain it is the workhorse loss method: it produces the excess rainfall that a unit hydrograph or reservoir routing scheme then transforms into a discharge hydrograph. This guide treats the method as an engineering procedure rather than a black box, covering the governing equations, the hydrologic soil group lookup, antecedent moisture adjustment, composite curve number assembly from spatial data, a production runoff function, validation, and the failure modes that trip up automated pipelines.

The method’s appeal is its parsimony. A single curve number, tabulated by the Natural Resources Conservation Service across decades of gauged small watersheds, captures the runoff response of a land parcel without requiring soil hydraulic conductivity, storage-depth profiles, or event-specific infiltration curves. That parsimony is also its limitation, and the failure-mode section below is explicit about where the abstraction breaks down. The two focused guides beneath this one handle the spatial side of the problem — computing composite curve numbers from land-use rasters in Python and estimating time of concentration with the NRCS velocity method — which together supply the two parameters a curve-number hydrograph model needs.

Prerequisites & Environment Setup

Curve number estimation is light on computation but heavy on spatial bookkeeping. The runoff equation is a two-line arithmetic expression; the effort goes into building a spatially consistent curve number grid from land cover and soils and then aggregating results to hydrologic units.

Input data requirements:

• A rainfall depth or a raster of rainfall depths, in consistent units (this guide uses millimetres)
• A curve number grid or a scalar CN, aligned to the rainfall grid when both are rasters
• Subbasin polygons for aggregation, in the same coordinate reference system as the rasters

Python stack:

conda create -n hydro-cn python=3.10 numpy rasterio geopandas rasterstats -c conda-forge
conda activate hydro-cn

Install the geospatial stack through conda-forge so that GDAL, PROJ, and GEOS resolve to a single compatible set of binaries. Mixed pip and conda installs are the most common source of silent CRS corruption when the CN grid and the rainfall grid are read through different GDAL builds. Because the runoff equation is applied cell by cell, both rasters must share an identical grid — same origin, resolution, and CRS — before the arithmetic runs.

Method Mechanics

The curve number method rests on a runoff-generation hypothesis: the ratio of actual retention to potential maximum retention equals the ratio of actual runoff to potential runoff (rainfall minus the initial abstraction). Rearranging that proportionality and substituting the initial abstraction gives the two equations at the core of the method.

Potential maximum retention. The curve number maps to a potential maximum retention S, the depth of water the soil column and surface can absorb once runoff begins. In millimetres:

S = 25400 / CN - 254

The equivalent imperial form, with S in inches, is S = 1000/CN - 10. A curve number of 100 gives S = 0 (every drop runs off, as for open water or pavement); a curve number of 50 gives S = 254 mm of retention. Because S is inversely and nonlinearly related to CN, small changes in CN at the high end (90 to 95) move runoff far more than the same change at the low end.

The runoff equation. With S known, the runoff depth Q for a storm of total depth P is:

Q = (P - 0.2S)^2 / (P + 0.8S)    when P > 0.2S
Q = 0                            when P <= 0.2S

The 0.2S term is the initial abstraction Ia — the rainfall intercepted, stored in surface depressions, and infiltrated before any runoff appears. The method fixes the initial abstraction ratio at 0.2, so Ia = 0.2S. A widely cited NRCS re-analysis of gauged events found that an initial abstraction ratio near 0.05 reproduces observed runoff better across many watersheds; adopting Ia = 0.05S requires rescaling the retention term and is discussed under failure modes. Unless you have recalibrated, keep the 0.2 ratio for consistency with published curve number tables and regulatory guidance from the NRCS.

The chart below shows the nonlinear rainfall-runoff response for three representative curve numbers. Note how a low CN both delays the onset of runoff (a larger initial abstraction) and depresses the runoff depth across the whole rainfall range.

Because the curve is convex, the marginal runoff coefficient — the fraction of an additional millimetre of rain that becomes runoff — rises with storm depth and approaches one for very large storms. This is physically reasonable: once the soil profile is saturated, nearly all further rainfall runs off. It is also why the method is more reliable for design storms than for small, frequent events, where the response is dominated by the initial abstraction term.

Hydrologic Soil Groups and the CN Lookup

The curve number for a parcel is read from a lookup table indexed by two attributes: the land cover or land use, and the hydrologic soil group (HSG) of the underlying soil. The HSG classifies soils A through D by their minimum infiltration rate when thoroughly wetted:

A given land cover generates progressively more runoff as the soil group moves from A to D. Row-crop agriculture on group A soil might carry a curve number near 67, while the same crop on group D soil is closer to 89. Open water and impervious surfaces sit at or near 98 regardless of the soil beneath them. Assembling the two-dimensional (land cover, HSG) -> CN table and applying it across a raster is the subject of the dedicated guide on computing composite curve numbers from land-use rasters in Python, which walks through pairing a land-cover raster with gridded soil survey data.

Antecedent Moisture Adjustment (AMC I / II / III)

The tabulated curve numbers assume average antecedent moisture — the amount of water already in the soil profile when the storm begins. This average state is antecedent moisture condition II (AMC II). Dry soils (AMC I) absorb more and produce less runoff for the same storm; saturated soils (AMC III) absorb less and produce more. Rather than maintain three separate lookup tables, the method converts the AMC II curve number to the drier or wetter condition with two closed-form relationships:

CN(I)   = 4.2 * CN(II) / (10 - 0.058 * CN(II))
CN(III) = 23  * CN(II) / (10 + 0.13  * CN(II))

For an AMC II curve number of 80, the dry-condition value drops to about 63 and the wet-condition value rises to about 91 — a swing that changes the runoff depth for a moderate storm by a factor of two or more. In continuous simulation these conversions are driven by a rolling five-day antecedent rainfall total; in single-event design work the condition is chosen to match the regulatory scenario, often AMC II for the standard case and AMC III for a conservative flood estimate. The conversions are monotonic and bounded, so they vectorize cleanly over a CN array.

Composite Curve Numbers

Real subbasins are heterogeneous: a single hydrologic unit spans several land covers over several soil groups. The standard practice collapses that mosaic into one representative curve number by area-weighting:

CN_composite = sum(area_i * CN_i) / sum(area_i)

where each i is a unique (land cover, HSG) polygon or raster cell within the subbasin. Weighting the curve numbers is the convention embedded in agency spreadsheets and most modeling software. It is exact only when the response is linear; because the runoff equation is convex, a composite CN slightly underestimates runoff for basins that mix very high and very low curve numbers under a large storm. For those cases, computing runoff per cell and then area-averaging the runoff depths — rather than the curve numbers — is the more defensible approach, and the production function below is written so that either aggregation order is available. The mechanics of zonal area-weighting with rasterstats are covered in the composite-CN guide.

Step-by-Step Workflow

A curve-number runoff estimate follows the same deterministic sequence whether it runs on a scalar CN or a continental raster.

Step 1 — Assemble the CN grid

Rasterize land cover and hydrologic soil group onto a common grid, then map each (land cover, HSG) pair to its tabulated curve number. The output is a CN raster at AMC II. Impervious fractions and open water must be handled explicitly so they are not diluted by an area-weighted average.

Step 2 — Apply the antecedent moisture adjustment

If the design scenario or the antecedent five-day rainfall calls for a condition other than AMC II, convert the whole grid with the closed-form relationships above. Keep the AMC II grid on disk so the adjustment is reproducible and auditable.

Step 3 — Convert CN to retention and compute runoff

Convert each cell’s curve number to retention S, compute the initial abstraction, and apply the runoff equation for the storm depth. The clamp on P - 0.2S is not optional: without it, cells where the storm is smaller than the initial abstraction produce a negative squared term that the equation would otherwise turn into spurious positive runoff.

Step 4 — Aggregate to subbasins

Reduce the cell-level runoff grid to a runoff depth per subbasin with a zonal mean, then hand those depths to a transform method. The excess this method produces is the direct input to the unit hydrograph methods that convert a depth of excess rainfall into a discharge hydrograph.

Production-Ready Code

The function below computes a runoff-depth grid from a single rainfall depth and a curve number array. It handles the AMC adjustment, the initial abstraction ratio as a parameter, nodata propagation, and the sub-threshold clamp. It is written to run on either a scalar CN or a full raster because both share the same vectorized arithmetic.

import logging
import numpy as np

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


def adjust_cn_for_amc(cn_ii: np.ndarray, amc: str) -> np.ndarray:
    """Convert an AMC II curve number array to AMC I (dry) or AMC III (wet)."""
    if amc == "II":
        return cn_ii
    if amc == "I":
        adjusted = 4.2 * cn_ii / (10.0 - 0.058 * cn_ii)
    elif amc == "III":
        adjusted = 23.0 * cn_ii / (10.0 + 0.13 * cn_ii)
    else:
        raise ValueError(f"amc must be 'I', 'II', or 'III'; got {amc!r}")
    logger.info("Applied AMC %s adjustment (mean CN %.1f -> %.1f)",
                amc, float(np.nanmean(cn_ii)), float(np.nanmean(adjusted)))
    return adjusted


def scs_runoff_depth(
    rainfall_mm: float,
    cn_array: np.ndarray,
    ia_ratio: float = 0.2,
    amc: str = "II",
    nodata: float = -9999.0,
) -> np.ndarray:
    """
    Compute SCS curve number runoff depth (mm) for one rainfall depth over a CN grid.

    Args:
        rainfall_mm: Storm rainfall depth in millimetres (scalar).
        cn_array:    Curve number array at AMC II, values in (0, 100]; nodata flagged.
        ia_ratio:    Initial abstraction ratio (0.2 standard; 0.05 requires rescaled S).
        amc:         Antecedent moisture condition 'I', 'II', or 'III'.
        nodata:      Sentinel for nodata cells, preserved in the output.

    Returns:
        Runoff depth array (mm), nodata cells set to the sentinel.
    """
    cn = np.asarray(cn_array, dtype="float64")
    valid = (cn > 0.0) & (cn <= 100.0)
    if not np.any(valid):
        raise ValueError("cn_array contains no valid curve numbers in (0, 100].")

    # AMC conversion on valid cells only, leaving nodata untouched
    cn_work = cn.copy()
    cn_work[valid] = adjust_cn_for_amc(cn[valid], amc)

    # Potential maximum retention S (mm); CN = 100 gives S = 0
    retention = np.full_like(cn, nodata)
    retention[valid] = 25400.0 / cn_work[valid] - 254.0

    initial_abstraction = ia_ratio * retention  # Ia = ratio * S

    # Runoff only where P exceeds Ia; clamp the numerator term at zero elsewhere
    runoff = np.full_like(cn, nodata)
    runs = valid & (rainfall_mm > initial_abstraction)
    excess = rainfall_mm - initial_abstraction[runs]
    denom = rainfall_mm + (1.0 - ia_ratio) * retention[runs]  # P + (1-ratio)*S
    runoff[runs] = excess ** 2 / denom

    # Valid cells that did not exceed Ia produce exactly zero runoff
    dry = valid & ~(rainfall_mm > initial_abstraction)
    runoff[dry] = 0.0

    logger.info(
        "Rainfall %.1f mm | AMC %s | Ia ratio %.2f | mean runoff %.2f mm over %d cells",
        rainfall_mm, amc, ia_ratio, float(np.mean(runoff[valid])), int(valid.sum())
    )
    return runoff


if __name__ == "__main__":
    # Demo: a 75 mm design storm over a small synthetic CN grid at AMC II
    demo_cn = np.array(
        [[70.0, 80.0, 90.0],
         [55.0, 85.0, 98.0],
         [-9999.0, 78.0, 62.0]],
        dtype="float64",
    )
    q_grid = scs_runoff_depth(rainfall_mm=75.0, cn_array=demo_cn, amc="II")
    logger.info("Runoff grid (mm):\n%s", np.round(q_grid, 1))

The denominator is written as P + (1 - ia_ratio) * S rather than the literal P + 0.8S so that the same function is correct when the initial abstraction ratio changes. With the default ratio of 0.2, (1 - 0.2) * S reduces to the familiar 0.8S.

Validation Protocol

Because the method is arithmetic, validation focuses on physical plausibility and spatial consistency rather than numerical convergence.

Runoff bounded by rainfall. For every cell, 0 <= Q <= P. A runoff depth exceeding the rainfall depth indicates a sign error or an unclamped sub-threshold term. Assert this bound explicitly after computing the grid.

Runoff coefficient sanity. The basin-average runoff coefficient Q/P should fall between about 0.1 for a permeable, low-CN watershed under a small storm and 0.9 for an impervious basin under a large storm. Coefficients outside that band point to a mis-assembled CN grid — often a land-cover class mapped to the wrong curve number.

Curve number distribution. Histogram the CN grid before computing runoff. Values below 30 or a large spike at 100 that does not correspond to mapped water or pavement usually mean the (land cover, HSG) lookup missed a class and fell back to a default. Cross-check the areal fractions against the land-cover source.

Independent event check. Where gauged runoff exists, compare the modeled runoff depth for observed storms against the measured direct runoff. Systematic over-prediction of small events is expected and is a property of the fixed 0.2 initial abstraction ratio, not a coding error; systematic bias across all event sizes indicates a CN calibration problem best addressed alongside the model calibration objective functions used to score the fit.

Common Failure Modes & Optimization

Sub-threshold storms and the P < Ia clamp. When the storm depth is smaller than the initial abstraction, P - 0.2S is negative and its square is positive — so an unclamped equation reports phantom runoff for storms that produce none. Always guard the equation with P > 0.2S and set runoff to zero otherwise. The production function does this by computing runoff only on the masked runs subset and assigning zero to the remaining valid cells.

Low-CN, low-rainfall instability. At low curve numbers the retention S is large, so the initial abstraction consumes most of a small storm and the denominator P + 0.8S is dominated by S. The runoff estimate in this regime is extremely sensitive to the CN value: a two-point CN error can double or halve the runoff. Treat low-CN results under small storms as order-of-magnitude estimates, and prefer per-cell runoff aggregation over composite-CN aggregation when a basin mixes low and high curve numbers.

Initial abstraction ratio choice. The 0.2 ratio is baked into the standard curve number tables. Switching to 0.05 without also rescaling S and the curve number produces internally inconsistent runoff, because the tabulated CN values were fit under the 0.2 assumption. If you adopt 0.05S, recalibrate the retention term against gauged events rather than reusing published curve numbers directly.

Grid misalignment. When the rainfall and CN grids differ in origin, resolution, or CRS, the cell-by-cell equation silently pairs mismatched locations. Resample both onto a shared grid first; the mechanics of reconciling projections are covered under coordinate reference system alignment. Verify identical shape and transform before the arithmetic runs.

Impervious dilution. Area-weighting a curve number over a subbasin that contains concentrated impervious area (a parking lot, a rooftop cluster) smears the high runoff of those surfaces across the whole unit and underestimates peak response. Model connected impervious area separately, or use a resolution fine enough that impervious cells retain their true curve number.

When to Use Curve Numbers vs. Alternatives

The curve number method is the right tool when you need an event-based runoff volume from readily available soils and land-cover data, when the watershed is ungauged, or when a regulatory framework mandates it. It is a loss method, not a full model: it tells you how much of a storm becomes runoff, not when that runoff arrives at the outlet.

To turn the excess rainfall into a hydrograph, pair the curve number loss with a transform. The unit hydrograph methods take the depth of excess this method produces and distribute it in time according to basin shape and response speed; the SCS dimensionless unit hydrograph in particular is the natural partner, since its timing is set by the same NRCS lag relationship. That timing depends on the basin’s time of concentration, computed in the companion guide on estimating time of concentration with the NRCS velocity method.

For situations where the curve number’s event-based abstraction is too coarse — long-duration continuous simulation, snowmelt-driven basins, or urban drainage with explicit conveyance — a physically based infiltration and routing engine is more appropriate. When you need to fit any of these models to observed discharge, the model calibration objective functions provide the Nash-Sutcliffe and Kling-Gupta scores used to judge whether a calibrated curve number reproduces reality.

Related Topics

• Computing composite curve numbers from land-use rasters in Python — building the (land cover, HSG) -> CN lookup and zonal area-weighting per subbasin
• Estimating time of concentration with the NRCS velocity method — deriving the Tc and lag that set the timing of the runoff this method produces
• Unit Hydrograph Methods — transforming the excess rainfall into a discharge hydrograph
• Model Calibration & Objective Functions — scoring a curve-number model against observed runoff
• Rainfall-Runoff Modeling & Hydrologic Simulation — the parent overview of loss methods, transforms, and calibration