In this study, we propose a regional Bayesian hierarchical model for flood frequency analysis. The Bayesian method is an alternative to the traditional regional flood frequency analysis. Instead of relying on the delineation of implicit homogeneous regions, the Bayesian hierarchical method describes the spatial dependence in its inner structure. Similar to the classical Bayesian hierarchical model, the process layer of our model presents the spatial variability of the parameters by considering different covariates (e.g., drainage area, elevation, precipitation). Beyond the three classical layers (data, process, and prior) of the Bayesian hierarchical model, we add a new layer referred to as the “L-moments layer”. The L-moments layer uses L-moments theory to select the best-fit probability distribution based on the available data. This new layer can overcome the subjective selection of the distribution based on extreme value theory and determine the distribution from the data instead. By adding this layer, we can combine the merits of regional flood frequency and Bayesian methods. A standard process of covariates selection is also proposed in the Bayesian hierarchical model. The performance of the Bayesian model is assessed by a case study over the Willamette River Basin in the Pacific Northwest, U.S. The uncertainty of different flood percentiles can be quantified from the posterior distributions using the Markov Chain Monte Carlo method. Temporal changes for the 100-year flood percentiles are also examined using a 20- and 30-year moving window method. The calculated shifts in flood risk can aid future water resources planning and management.