Estimating spatially distributed hydrological model parameters on ungauged catchments poses a challenging regionalization problem, especially when searching for a transfer function that relates quantitatively physical descriptors to conceptual hydrological parameters, and imposing spatial constraints needed given sparse constraining discharge data. This paper introduces a Hybrid Variational Data Assimilation Parameter Regionalization (HVDA-PR) approach extending the multiscale parameter regionalization (MPR) technique. HVDA-PR leverages spatially distributed cost gradients to infer complex transfer functions designed for high-resolution hydrological models. The key components of HVDA-PR involve incorporating learnable regionalization mappings, which consist of either multivariate regressions or neural networks, into a differentiable hydrological model. This enables the exploitation of the informative content of heterogeneous datasets across extensive spatio-temporal computational domains, particularly in high-dimensional regionalization, with adapted optimization algorithms and accurate adjoint-based gradients. The inverse problem was tackled with a multi-gauge calibration cost function accounting for information from multiple observation sites. HVDA-PR was tested on high-resolution, hourly and kilometric regional modeling of two flash flood prone study areas located in the South of France. In both study areas, the median NSE scores of HVDA-PR ranged from 0.52 to 0.78 at pseudo-ungauged sites over calibration and validation periods, which exhibited strong regionalization performance, improving NSE by up to 0.57 compared to the baseline regionalization model calibrated with lumped parameters, and achieving a comparable level to the reference solution obtained from local uniform calibration (median NSE from 0.59 to 0.79). Multiple validation metrics based on hydrological flood-scale signatures are employed to assess the accuracy and robustness of the approach. The ability to produce physically explainable parameter maps from the physical descriptors demonstrates the effectiveness of HVDA-PR, in addition to the impressive scores obtained in calibration, and in validation in time and at ungauged sites. The regionalization method is amenable to state-parameter correction from multi-source data, at multiple time scales such as for operational data assimilation, and it is transposable to other differentiable geophysical models.