Computing Nash-Sutcliffe Efficiency and KGE in Python

Nash-Sutcliffe efficiency (NSE) and Kling-Gupta efficiency (KGE) are the two scalar scores that decide whether a calibrated hydrograph is acceptable, and getting their implementation exactly right matters because both are easy to compute almost correctly. This page is the metric reference for the Model Calibration & Objective Functions workflow, itself part of the broader Rainfall-Runoff Modeling & Hydrologic Simulation domain. It gives the exact formulas, numpy/pandas implementations that return every metric plus the KGE component decomposition, the published interpretation thresholds, and the alignment and NaN handling that silently corrupt these numbers when skipped.

Everything below assumes two discharge series — one simulated, one observed from a gauge such as USGS NWIS — sampled at the same timestep and covering the same period.

The Formulas

Nash-Sutcliffe Efficiency

NSE normalizes the model’s squared error by the variance of the observations:

NSE = 1 − [ Σ (Q_obs,t − Q_sim,t)² ] / [ Σ (Q_obs,t − mean(Q_obs))² ]

The denominator is the total sum of squares of the observations — their variance times the sample count. NSE = 1 is a perfect match. NSE = 0 means the model’s squared error equals the observed variance, i.e. the observed mean is exactly as good a predictor. NSE < 0 means the mean is a better predictor. Because errors are squared, a handful of large peaks dominate the statistic and low-flow error is nearly invisible.

Kling-Gupta Efficiency (Gupta et al., 2009)

KGE measures the Euclidean distance from the ideal point (1, 1, 1) in the three-dimensional space of correlation, variability ratio, and bias ratio:

KGE = 1 − sqrt( (r − 1)² + (α − 1)² + (β − 1)² )

r = Pearson correlation( Q_sim , Q_obs )
α = std(Q_sim) / std(Q_obs)      # variability ratio
β = mean(Q_sim) / mean(Q_obs)    # bias ratio

KGE = 1 is perfect. The three terms are independent: a model can have perfect correlation (r = 1) yet still be penalized for under-dispersion (α < 1) or volume bias (β ≠ 1). This separation is exactly why KGE is preferred over NSE as a single calibration criterion — NSE conflates these errors and its optimum systematically damps variability.

PBIAS

Percent bias measures the tendency of the simulated volumes to be larger or smaller than observed:

PBIAS = 100 × [ Σ (Q_obs,t − Q_sim,t) ] / [ Σ Q_obs,t ]

With this (Moriasi et al., 2007) sign convention, positive PBIAS = model under-prediction of total volume and negative PBIAS = over-prediction. PBIAS = 0 is ideal. Some libraries flip the sign by using sim − obs in the numerator, so always state the convention when you report it.

logNSE

logNSE is ordinary NSE computed on log-transformed flows. The log transform compresses the peaks and stretches the low-flow range, moving the metric’s sensitivity from floods to baseflow and recession:

logNSE = NSE( ln(Q_obs + ε) , ln(Q_sim + ε) )

The offset ε prevents ln(0) on zero-flow days. A standard choice is ε = 0.01 × mean(Q_obs) — small enough not to distort the low flows it is meant to preserve, large enough to keep the logarithm finite.

Interpretation Thresholds

Published performance ratings (Moriasi et al., 2007, for monthly streamflow, with KGE benchmarks after Knoben et al., 2019) give defensible acceptance cutoffs. Sub-daily or event-scale work is typically judged more leniently.

Two benchmark values are easy to confuse. For NSE, 0 is the mean-flow benchmark — a model at NSE = 0 is no better than the observed mean. For KGE, the mean-flow benchmark is approximately −0.41, so any positive KGE already beats the mean. Judging a KGE model against a threshold of 0 as if it were NSE unfairly rejects usable models.

Annotated Implementation

The function below returns every metric in one pass and exposes the KGE components so a poor score can be diagnosed rather than merely reported.

import logging
import numpy as np
import pandas as pd

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


def streamflow_metrics(
    observed: pd.Series,
    simulated: pd.Series,
    log_epsilon_frac: float = 0.01,
) -> dict:
    """
    Compute NSE, KGE (with r/alpha/beta), PBIAS, and logNSE for two aligned
    discharge series.

    Parameters
    ----------
    observed, simulated : pandas.Series
        Discharge with a shared DatetimeIndex. Series are inner-joined on their
        index and timesteps missing in either series are dropped.
    log_epsilon_frac : float
        Fraction of mean observed flow added before the log transform, so that
        zero-flow days do not produce -inf in logNSE.

    Returns
    -------
    dict of metric name -> float
    """
    # --- 1. Align on the shared index and drop NaNs from EITHER series ---
    paired = pd.concat({"obs": observed, "sim": simulated}, axis=1).dropna()
    n = len(paired)
    if n < 2:
        raise ValueError(f"Need >=2 overlapping timesteps; got {n}.")
    obs = paired["obs"].to_numpy(dtype=float)
    sim = paired["sim"].to_numpy(dtype=float)
    logger.info("Scoring %d aligned timesteps", n)

    # --- 2. NSE ---
    obs_variance = np.sum((obs - obs.mean()) ** 2)
    if obs_variance == 0:
        raise ValueError("Observed series has zero variance; NSE/KGE undefined.")
    nse = 1.0 - np.sum((obs - sim) ** 2) / obs_variance

    # --- 3. KGE and its three components ---
    # Pearson r; guard against a constant simulated series (std == 0)
    if sim.std(ddof=0) == 0:
        r = 0.0
        logger.warning("Simulated series is constant; correlation set to 0.")
    else:
        r = np.corrcoef(sim, obs)[0, 1]
    alpha = sim.std(ddof=0) / obs.std(ddof=0)      # variability ratio
    beta = sim.mean() / obs.mean()                 # bias ratio
    kge = 1.0 - np.sqrt((r - 1) ** 2 + (alpha - 1) ** 2 + (beta - 1) ** 2)

    # --- 4. PBIAS (Moriasi sign: positive => under-prediction of volume) ---
    pbias = 100.0 * np.sum(obs - sim) / np.sum(obs)

    # --- 5. logNSE with an epsilon offset ---
    eps = log_epsilon_frac * obs.mean()
    log_obs = np.log(obs + eps)
    log_sim = np.log(sim + eps)
    log_denom = np.sum((log_obs - log_obs.mean()) ** 2)
    log_nse = 1.0 - np.sum((log_obs - log_sim) ** 2) / log_denom

    metrics = {
        "NSE": float(nse),
        "KGE": float(kge),
        "KGE_r": float(r),
        "KGE_alpha": float(alpha),
        "KGE_beta": float(beta),
        "PBIAS_pct": float(pbias),
        "logNSE": float(log_nse),
        "n": n,
    }
    logger.info(
        "NSE=%.3f KGE=%.3f (r=%.3f a=%.3f b=%.3f) PBIAS=%.1f%% logNSE=%.3f",
        metrics["NSE"], metrics["KGE"], metrics["KGE_r"],
        metrics["KGE_alpha"], metrics["KGE_beta"],
        metrics["PBIAS_pct"], metrics["logNSE"],
    )
    return metrics


if __name__ == "__main__":
    rng = np.random.default_rng(7)
    idx = pd.date_range("2016-01-01", periods=730, freq="D")
    obs = pd.Series(np.abs(rng.gamma(1.5, 3.0, size=730)) + 0.2, index=idx)
    # A simulation that lags and slightly under-disperses the observations
    sim = pd.Series(0.9 * obs.shift(1).fillna(obs.iloc[0]).to_numpy() + 0.3, index=idx)
    result = streamflow_metrics(obs, sim)
    print(result)

Reading the KGE decomposition

The returned KGE_r, KGE_alpha, and KGE_beta turn the aggregate score into a diagnosis:

• r below 1 — timing or shape error. The rising limbs and peaks do not line up. Look at routing lag, the unit hydrograph, or timestep offsets.
• alpha below 1 — the simulation is too smooth; its flow variability is smaller than observed. Peaks are clipped and recessions too flat.
• alpha above 1 — the simulation is too flashy.
• beta away from 1 — volume bias. beta = 1.10 means the model produces 10% too much total runoff; this maps directly onto a PBIAS problem and usually a water-balance parameter (loss rate, curve number, evapotranspiration).

A KGE of 0.6 with beta = 1.0 and alpha = 0.7 is a variability problem, not a bias problem — a distinction the single number cannot convey and the reason the automated search in automated hydrologic model calibration with SciPy and SPOTPY benefits from a decomposed objective.

A worked interpretation

Suppose streamflow_metrics returns NSE = 0.71, KGE = 0.68 with r = 0.94, alpha = 0.78, beta = 1.02, PBIAS = -1.9%, and logNSE = 0.41. Read left to right. The NSE of 0.71 and KGE of 0.68 both fall in the “good” band, and beta = 1.02 with PBIAS near −2% confirms the water balance is essentially correct — the model produces about 2% too much total volume. The r of 0.94 says timing and shape are strong. The story is entirely in alpha = 0.78: the simulation carries only 78% of the observed flow variability, so its peaks are clipped and its recessions are too flat. Finally, the logNSE of 0.41 — much lower than the linear NSE — reveals that the low-flow fit is only marginal even though the peak-dominated NSE looks healthy. The correct next move is not a volume adjustment (the balance is fine) but a change that sharpens peaks and steepens recessions, and a re-check of the baseflow parameters that logNSE is flagging. Without the decomposition, an NSE of 0.71 alone would have suggested the model was simply “good” and hidden both the under-dispersion and the weak low-flow behaviour.

Gotchas That Corrupt These Metrics

Timestep misalignment. If observed is daily and simulated is hourly, a naive element-wise subtraction compares unrelated timesteps and produces a meaningless score. Resample both to a common frequency (simulated.resample("D").mean()) before pairing, and confirm the two indices are identical afterward. The pd.concat(...).dropna() pattern above aligns on the index rather than on position, which is the safe default.

NaN handling. A single NaN in either series poisons np.sum and np.corrcoef, returning NaN for the whole metric. Gauge records routinely contain gaps (ice-affected periods, sensor outages), so always drop timesteps missing in either series — not just the simulated one — before computing anything. Dropping from only one side leaves length-mismatched arrays.

Zero observed variance. If the observed window is flat (a fully regulated reach, or a too-short dry-season slice), the NSE and KGE denominators are zero and both metrics are undefined. The implementation raises rather than returning a silent inf; catch this and lengthen or move the evaluation window.

Why the KGE decomposition matters. Reporting only the scalar KGE discards the information needed to fix the model. Two calibrations with identical KGE = 0.65 can require opposite corrections — one has a volume bias (beta off), the other a variability problem (alpha off). Always carry r, alpha, and beta through to the report.

Log offset choice. Setting ε too large in logNSE flattens the low-flow signal the metric exists to capture; setting it to zero throws on any zero-flow day. Tie it to mean flow (a small percentage) rather than using a fixed absolute constant, so the offset scales with the basin.

Related Topics

• Model Calibration & Objective Functions — the parent workflow that uses these metrics as the loss the optimizer minimizes
• Automated Hydrologic Model Calibration with SciPy and SPOTPY — the sibling guide that feeds these objective functions into a global optimizer
• Rainfall-Runoff Modeling & Hydrologic Simulation — the domain overview covering model building, simulation, and calibration end to end