Automated Hydrologic Model Calibration with SciPy and SPOTPY

Automatic calibration replaces manual parameter tweaking with a global optimizer that proposes parameter sets, runs the model, and follows the objective downhill until it converges. This guide implements that search two ways — with scipy.optimize.differential_evolution and with the SCE-UA algorithm through a spotpy setup class — as a focused how-to under the Model Calibration & Objective Functions topic, which sits inside the Rainfall-Runoff Modeling & Hydrologic Simulation domain. Both optimizers minimize the same loss; the objective they minimize is one of the efficiency metrics defined in computing Nash-Sutcliffe efficiency and KGE in Python.

The three pieces every automatic calibration needs are the same regardless of optimizer: a set of parameter bounds, an objective wrapper that turns a parameter vector into a loss, and a search algorithm. This page builds each, for both tools.

Prerequisites

You need a model exposed as a callable run_model(params) -> simulated_series, an observed discharge series (for example a USGS NWIS gauge), and physically defensible bounds for each free parameter. Install the two optimizers:

pip install scipy spotpy

differential_evolution ships with SciPy and needs no extra dependency. spotpy (Statistical Parameter Optimization Tool for Python) adds SCE-UA, DREAM, Latin-hypercube sampling, and the analyser utilities that read a results database.

Step 1 — Parameter Bounds and the Objective Wrapper

The objective wrapper is the contract shared by every optimizer: it accepts a parameter vector, runs the model, scores the result, and returns a scalar loss to minimize. Because both NSE and KGE are maximized at 1, the wrapper returns 1 − metric.

import logging
import numpy as np
import pandas as pd

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

# (low, high) bounds in parameter order: e.g. loss rate, recession k, routing lag
PARAM_BOUNDS = [(0.0, 50.0), (0.50, 0.99), (0.0, 24.0)]
PARAM_NAMES = ["loss_mm", "recession_k", "lag_hr"]


def objective_loss(params, run_model, observed, cal_slice) -> float:
    """Run the model for `params` and return 1 - KGE over the calibration window."""
    simulated = run_model(np.asarray(params))
    paired = pd.concat(
        {"obs": observed.loc[cal_slice], "sim": simulated.loc[cal_slice]}, axis=1
    ).dropna()
    if len(paired) < 2:
        logger.warning("Insufficient overlap for params=%s", np.round(params, 3))
        return 1e6
    obs, sim = paired["obs"].to_numpy(), paired["sim"].to_numpy()
    r = np.corrcoef(sim, obs)[0, 1] if sim.std() > 0 else 0.0
    alpha = sim.std(ddof=0) / obs.std(ddof=0)
    beta = sim.mean() / obs.mean()
    kge = 1.0 - np.sqrt((r - 1) ** 2 + (alpha - 1) ** 2 + (beta - 1) ** 2)
    return 1.0 - kge

Bounds do double duty: they constrain the search to physically meaningful values and they define the dimensionality the optimizer explores. Every parameter you add roughly multiplies the search cost, so include only the ones a sensitivity screen flagged as identifiable.

Step 2 — SciPy differential_evolution

differential_evolution is a population-based global optimizer. It maintains a population of popsize × n_params candidate vectors, mutates and recombines them each generation, and keeps whichever improves the loss. It needs no gradient, tolerates the discontinuous loss surfaces of hydrologic models, and lives entirely in SciPy.

from functools import partial
from scipy.optimize import differential_evolution


def calibrate_with_de(run_model, observed, cal_slice,
                      popsize=15, maxiter=300, tol=1e-3, seed=42):
    """Global calibration with SciPy differential evolution."""
    loss = partial(objective_loss, run_model=run_model,
                   observed=observed, cal_slice=cal_slice)

    logger.info("differential_evolution: %d params, popsize=%d, maxiter=%d",
                len(PARAM_BOUNDS), popsize, maxiter)
    result = differential_evolution(
        loss,
        bounds=PARAM_BOUNDS,
        popsize=popsize,
        maxiter=maxiter,
        tol=tol,
        seed=seed,            # reproducibility
        mutation=(0.5, 1.0),  # dithered mutation factor
        recombination=0.7,
        polish=True,          # local refinement of the best vector at the end
        updating="deferred",  # evaluate a whole generation at once
        workers=-1,           # use all CPU cores (needs deferred updating)
    )
    best = dict(zip(PARAM_NAMES, result.x))
    logger.info("Best KGE=%.4f at %s (%d evals)",
                1.0 - result.fun, {k: round(v, 3) for k, v in best.items()},
                result.nfev)
    return result.x, 1.0 - result.fun

result.x is the best parameter vector and result.fun is its loss, so 1 − result.fun is the calibrated KGE. result.nfev reports how many model runs the search consumed — the number to watch when each run is expensive.

Step 3 — SPOTPY SCE-UA

SCE-UA (Shuffled Complex Evolution, University of Arizona) was designed specifically for conceptual rainfall-runoff calibration and is the field’s reference global optimizer. spotpy drives it through a setup class that implements four methods: parameters, simulation, evaluation, and objectivefunction.

import spotpy
from spotpy.parameter import Uniform


class HydroSpotSetup:
    """spotpy setup class wrapping a rainfall-runoff model for SCE-UA calibration."""

    def __init__(self, run_model, observed, cal_slice):
        self.run_model = run_model
        self.cal_slice = cal_slice
        self.observed = observed.loc[cal_slice]
        # One Uniform parameter per free variable, with bounds from PARAM_BOUNDS
        self._params = [
            Uniform(name, low=lo, high=hi)
            for name, (lo, hi) in zip(PARAM_NAMES, PARAM_BOUNDS)
        ]

    def parameters(self):
        # spotpy calls this to draw a candidate vector each iteration
        return spotpy.parameter.generate(self._params)

    def simulation(self, vector):
        # Run the model for the sampled parameter vector; return simulated flow
        sim = self.run_model(np.asarray(list(vector)))
        return sim.loc[self.cal_slice].to_numpy()

    def evaluation(self):
        # The observations the simulation is scored against
        return self.observed.to_numpy()

    def objectivefunction(self, simulation, evaluation):
        # spotpy MINIMIZES by default for sceua; return 1 - KGE
        kge = spotpy.objectivefunctions.kge(evaluation, simulation)
        return 1.0 - kge


def calibrate_with_sceua(run_model, observed, cal_slice,
                         repetitions=8000, ngs=None, dbname="sceua_hydro"):
    """Global calibration with SPOTPY's SCE-UA sampler."""
    setup = HydroSpotSetup(run_model, observed, cal_slice)
    ngs = ngs or (2 * len(PARAM_BOUNDS) + 1)   # complexes; a common default
    logger.info("SCE-UA: repetitions<=%d, complexes(ngs)=%d", repetitions, ngs)

    sampler = spotpy.algorithms.sceua(
        setup, dbname=dbname, dbformat="csv", parallel="seq",
    )
    sampler.sample(repetitions, ngs=ngs)

    # --- Read the best parameter set from the results database ---
    results = spotpy.analyser.load_csv_results(dbname)
    best_index = spotpy.analyser.get_minlikeindex(results)[0]  # min loss row
    best_row = results[best_index]
    best_params = np.array([best_row[f"par{n}"] for n in PARAM_NAMES])
    best_kge = 1.0 - best_row["like1"]
    logger.info("SCE-UA best KGE=%.4f at %s",
                best_kge, dict(zip(PARAM_NAMES, np.round(best_params, 3))))
    return best_params, best_kge

get_minlikeindex returns the row with the smallest loss because the setup minimizes 1 − KGE; spotpy.analyser.get_best_parameterset is the higher-level convenience that does the same lookup. The results CSV also holds every sampled set, which is what you feed into an uncertainty analysis when you want more than a point estimate.

Optimizer Settings Reference

A practical starting point for a 3–6 parameter conceptual model is popsize=15, maxiter=200–300 for differential evolution, or repetitions=5000–15000 with ngs = 2·n_params + 1 complexes for SCE-UA. Increase both if restarts from different seeds land on different optima.

The two population controls trade coverage against cost in the same way. A larger popsize or ngs samples more of the parameter space per generation, which lowers the risk of the search converging into a local minimum but multiplies the number of model runs. The iteration controls (maxiter, repetitions) bound how long the refinement continues; because SCE-UA also stops on its internal convergence tests (peps, the parameter-convergence tolerance, and pcento, the objective-change tolerance), its repetitions value behaves as a ceiling rather than a fixed run count. Differential evolution has no early-stop of that kind beyond tol, so its maxiter × popsize × n_params product is closer to the actual worst-case run budget you should plan for.

Parallel Evaluation

Because every candidate parameter set requires an independent model run, calibration is embarrassingly parallel within a generation. SciPy exposes this through workers=-1, which forks the population evaluation across all available cores; it requires updating="deferred" so that the whole generation is scored before the next is proposed. spotpy distributes runs with parallel="mpc" (multiprocessing on one machine) or parallel="mpi" (across nodes under an MPI runtime such as mpirun). In both cases the speed-up is bounded by the slowest run in each generation and by the fixed cost of forking or serializing the model wrapper, so a 16-core machine rarely delivers a full 16× speed-up — 8–10× is a realistic expectation for second-scale model runs. The hard prerequisite is a picklable, stateless wrapper: if run_model writes intermediate files, give each evaluation a unique temporary directory so concurrent workers never collide on the same path.

Reading and Validating the Result

Whichever optimizer you run, the deliverable is a best parameter vector plus its calibration score. Do not stop there:

1. Re-run the model at the best vector over the full record and compute the full metric set (NSE, KGE with its components, PBIAS, logNSE) from the NSE and KGE reference.
2. Score on a holdout window the optimizer never saw, and report the calibration-to-validation drop.
3. Restart from independent seeds. If two seeds converge to materially different parameter vectors with similar scores, you are seeing equifinality and should carry an ensemble rather than a single vector.

Gotchas

Expensive evaluations dominate wall-clock time. Each iteration runs the full model, so a search of 10,000 evaluations at 2 seconds per run is more than five hours single-threaded. Profile one model run first, then size maxiter/repetitions against your time budget, and parallelize before scaling up.

Equifinality produces unstable optima. Different seeds landing on different “best” parameters with near-identical scores is the normal signature of an over-parameterized model, not an optimizer bug. Screen out insensitive parameters and, where uncertainty matters, switch from a single-best sampler to spotpy’s DREAM to obtain a posterior distribution.

Unseeded runs are irreproducible. Both optimizers are stochastic. Always fix seed (differential evolution) or the NumPy RNG (spotpy) so a reviewer can reproduce the calibration exactly.

Bounds that are too wide or too tight. Over-wide bounds waste evaluations in physically absurd regions and worsen equifinality; over-tight bounds can exclude the true optimum. Derive bounds from physical reasoning and widen only if the optimum pins against an edge.

Parallelism needs picklable, stateless models. workers=-1 and parallel="mpc"/"mpi" fork the objective across processes. If run_model writes to a shared file or holds global state, concurrent evaluations corrupt each other. Make the wrapper side-effect-free — write to per-run temporary paths — before enabling parallel evaluation.

Related Topics

• Model Calibration & Objective Functions — the parent workflow covering manual vs automatic calibration, validation, and equifinality
• Computing Nash-Sutcliffe Efficiency and KGE in Python — the objective functions these optimizers minimize, with formulas and thresholds
• Rainfall-Runoff Modeling & Hydrologic Simulation — the domain overview linking model building, simulation, and calibration