Recently, runoff simulations in small, headwater basins have been improved by methodological advances such as deep learning (DL). Hydrologic routing modules are typically needed to simulate flows in stem rivers downstream of large, heterogeneous basins, but obtaining suitable parameterization for them has previously been difficult. It is unclear if downstream daily discharge contains enough information to constrain spatially-distributed parameterization. Building on recent advances in differentiable modeling principles, here we propose a differentiable, learnable physics-based routing model. It mimics the classical Muskingum-Cunge routing model but embeds a neural network (NN) to provide parameterizations for Manning’s roughness coefficient (n) and channel geometries. The embedded NN, which uses (imperfect) DL-simulated runoffs as the forcing data and reach-scale attributes as inputs, was trained solely on downstream hydrographs. Our synthetic experiments show that while channel geometries cannot be identified, we can learn a parameterization scheme for n that captures the overall spatial pattern. Training on short real-world data showed that we could obtain highly accurate routing results for both the training and inner, untrained gages. For larger basins, our results are better than a DL model assuming homogeneity or the sum of runoff from subbasins. The parameterization learned from a short training period gave high performance in other periods, despite significant bias in runoff. This is the first time an interpretable, physics-based model is learned on the river network to infer spatially-distributed parameters. The trained n parameterization can be coupled to traditional runoff models and ported to traditional programming environments.