This study presents a method to improve the reconstruction of historical flows on gauged and ungauged basins. To do so, a multi-model weighted averaging of hydrological model simulations for a large spatial domain in the province of Quebec, Canada, is used. The distributed hydrological model HYDROTEL was implemented over the region and covered 95 gauged basins. An optimal interpolation (OI) assimilation method was first implemented as a baseline to improve the HYDROTEL flow simulations over the 95 basins. Then, a series of multi-model averaging techniques were applied to an ensemble of 144 HYDROTEL simulations that were generated by modifying parameter sets, driving weather datasets and evapotranspiration modules. The averaging methods were applied in a leave-one-out cross-validation scheme, where all 94 gauged basins were pooled together to compute the weights, and those weights were applied to the 95th basin. All basins were evaluated in such a manner and compared to the OI method. Implementing a year-by-year (or shorter period) weighting scheme instead of computing weights over all available data significantly improved the results. This allowed the weights to better reflect each year’s hydrological characteristics rather than compromising to improve the overall average. The Kling-Gupta Efficiency (KGE) and peak flow metrics showed that the Granger-Ramanathan variant “A” (GRA) was similar in performance to the OI method but did not have the drawbacks that OI can typically introduce. The multi-model application can also be further improved by adding more simulations from other hydrological models, whereas the OI method cannot make use of such additional information, thus hitting a performance plateau. This study shows that it is possible to improve regional hydrological model simulations, both for overall flows and peak flows, and on the historical period for both gauged and ungauged basins. This can then be used to better estimate risk in flood frequency analysis and other statistical analyses.