Flood frequency analysis at large scales, essential for the development of flood risk maps, is hindered by the uncertainty in outputs of process-based hydrological models and the scarcity of gauge flow data. We develop a Bayesian hierarchical model (BHM) based on the Generalized Extreme Value (GEV) distribution for regional flood frequency analysis at high resolution across North America. Our model leverages annual maximum flow data from ≈20000 gauged stations and a dataset of 130 static catchment-specific covariates to predict extreme flows at all catchments over the continent as well as their associated statistical uncertainty. Additionally, a modification is made to the data layer of the BHM to include daily gauge discharge data when available, which improves the precision of the discharge level estimates. We validated the model using a hold-out approach and found that its predictive power is very good for the GEV distribution location and scale parameters and improveable for the shape parameter, which is notoriously hard to estimate. The resulting discharge return levels yield satisfying agreement when compared with available design peak discharge from various government sources. Assessment of the covariates’ contributions to the model also informs on the most relevant underlying factors influencing flood-inducing peak flows. The key covariates in our model are temperature-related bioindicators, the catchment drainage area and geographical location.