Reclaiming First Principles: A Differentiable Framework for Conceptual Hydrologic Models

TL;DR

This paper develops an exact analytic framework for differentiable hydrologic modeling by augmenting the governing ODEs with forward sensitivity equations to produce state-and-parameter Jacobians entirely analytically. The resulting gradients, computed as $\mathbf{g}_{n}(\boldsymbol{\uptheta}) = \mathbf{J}_{q}(\boldsymbol{\uptheta})^{\top} \boldsymbol{\updelta}_{n}(\boldsymbol{\uptheta})$, enable gradient-based calibration across diverse differentiable losses (e.g., $L_{1}$, $L_{2}$, NSE, KGE, Huber, FDC) without finite-difference noise or autodiff overhead. The method is demonstrated on Hymod, HMODEL, SACSMA, and Xinanjiang, showing near-perfect agreement with numerical Jacobians while delivering dramatic CPU-time speedups (roughly $70$–$500\times$) over finite differences and orders of magnitude faster gradients than automatic differentiation for long time series. The results jointly offer faster, more stable calibration, transparent sensitivity analyses, and insights into parameter identifiability and equifinality, with clear pathways to large-scale, data-driven hydrologic learning and uncertainty quantification. The framework preserves physical interpretability and can be embedded directly in existing hydrologic codes, avoiding opaque or CPU-intensive autodiff toolchains.

Abstract

Conceptual hydrologic models remain the cornerstone of rainfall-runoff modeling, yet their calibration is often slow and numerically fragile. Most gradient-based parameter estimation methods rely on finite-difference approximations or automatic differentiation frameworks (e.g., JAX, PyTorch and TensorFlow), which are computationally demanding and introduce truncation errors, solver instabilities, and substantial overhead. These limitations are particularly acute for the ODE systems of conceptual watershed models. Here we introduce a fully analytic and computationally efficient framework for differentiable hydrologic modeling based on exact parameter sensitivities. By augmenting the governing ODE system with sensitivity equations, we jointly evolve the model states and the Jacobian matrix with respect to all parameters. This Jacobian then provides fully analytic gradient vectors for any differentiable loss function. These include classical objective functions such as the sum of absolute and squared residuals, widely used hydrologic performance metrics such as the Nash-Sutcliffe and Kling-Gupta efficiencies, robust loss functions that down-weight extreme events, and hydrograph-based functionals such as flow-duration and recession curves. The analytic sensitivities eliminate the step-size dependence and noise inherent to numerical differentiation, while avoiding the instability of adjoint methods and the overhead of modern machine-learning autodiff toolchains. The resulting gradients are deterministic, physically interpretable, and straightforward to embed in gradient-based optimizers. Overall, this work enables rapid, stable, and transparent gradient-based calibration of conceptual hydrologic models, unlocking the full potential of differentiable modeling without reliance on external, opaque, or CPU-intensive automatic-differentiation libraries.

Figures (13)

• Figure 1: Schematic illustration of the parameter transformations used for differentiable hydrologic modeling. The left panel depicts a watershed together with its original, physically interpretable hydrologic parameters $\boldsymbol{\uptheta} = (\theta_{1},\ldots,\theta_{d})^{\top}$. These parameters are first mapped to the unit hypercube, producing $\underline{\boldsymbol{\uptheta}} = (\underline{\theta}_{1},\ldots,\underline{\theta}_{d})^{\top} \in (0,1)^{d}$, shown in the center panel. The right panel displays the corresponding unconstrained parameter vector $\boldsymbol{\upvartheta} = (\vartheta_{1},\ldots,\vartheta_{d})^{\top} \in \mathbb{R}^{d}$, obtained via the logit transformation $\underline{\theta}_{j} = 1/\{(1 + \exp(-\vartheta_{j})\}$. The inverse mapping, $\vartheta_{j} = \log(\underline{\theta}_{j}/\{1 -\underline{\theta}_{j}\})$, returns parameters from the unconstrained space to the unit cube $\mathbb{U}^{d}$. This reparameterization eliminates boundary constraints yet maintains a smooth, invertible mapping to the physical hydrologic parameter space. As a result, gradient-based calibration becomes both more efficient and more robust.
• Figure 2: Schematic illustration of Richardson extrapolation for numerical differentiation. Starting from an anchor point $x_{0}$ (blue square), the function $f(x)$ (gray curve) is evaluated at a sequence of logarithmically spaced perturbations $x_{0} \pm h_{j}$ (black markers) using a central differencing scheme. Each perturbation pair yields a finite-difference approximation of the derivative, $f'(x_{0}) \approx \frac{1}{2}h^{-1}_{j}\{f(x_{0} + h_{j}) - f(x_{0}-h_{j})\}$ at $x_{0}$, illustrated with the colored secant lines for step sizes $h_{1}$ and $h_{5}$. The derivative estimates for the progressively smaller step sizes $h_{j}$, are recursively combined using multi-term Richardson (Romberg) extrapolation to cancel leading-order truncation errors. This procedure provides a high-order approximation of the derivative in the limit $h \rightarrow 0$, indicated by the brown line.
• Figure 3: Sensitivity of simulated streamflow $q_{1}, \ldots, q_{n}$ to xinanjiang parameters (a) $f_\mathrm{p}$, (b) $A_\mathrm{im}$, (c) $a$, (d) $b$, (e) $f_\mathrm{wm}$, (f) $f_\mathrm{lm}$, (g) $c$, (h) $s_\mathrm{tot}$, (i) $\beta$, (j) $k_\mathrm{i}$, (k) $k_\mathrm{g}$, and (l) $c_\mathrm{i}$ for the Wye river, UK. The entries of the different parameters make up the Jacobian matrix $\mathbf{J}_{q}(\boldsymbol{\uptheta}) = \partial \mathbf{q}_{n}/\partial \boldsymbol{\uptheta}^{\top}$.
• Figure 4: Comparison of analytic and numerical Jacobian entries for (a) hymod ( $786,640$ entries), (b) hmodel ( $1,101,296$ entries), (c) sacsma ($2,045,264$ entries), and (d) xinanjiang ($2,202,592$ entries) based on $N = 10$ randomly sampled parameter vectors across the Leaf River, French Broad, Wye, and Severn catchments. In all cases, points fall tightly along the 1:1 line, demonstrating near-perfect agreement between analytic sensitivities and their numerical finite-difference counterparts computed using the DERIVESTsuite toolbox of derrico2024.
• Figure 5: Comparison of analytic and numerical gradient vectors for four conceptual hydrologic models using $N = 10$ randomly sampled parameter vectors across the Leaf River, French Broad, Wye, and Severn catchments. Each panel shows a scatter plot of numerical versus analytic gradient values for (a) hymod ($960$ entries), (b) hmodel ($1,344$ entries), (c) sacsma ($2,496$ entries), and (d) xinanjiang ($2,688$ entries).