As atmospheric moisture capacity is highly sensitive to rising temperatures, precipitation extremes are widely projected to intensify with a warming climate and thus altering the flooding generation regime. Previous works seldomly focused on bivariate flood quantiles under climate change at a national scale, and fewer flooding projections quantified the estimation uncertainty sourced from sample size limitation. This study systematically investigates the changes in bivariate quantiles of flood peak and volume with incorporation of sampling uncertainty for 151 catchments over China, with climate trajectories projected by a set of multi-model ensemble under representative concentration pathway (RCP) 8.5. After correcting the systematical biases of eight CMIP5 GCM outputs, four state-of-the-art hydrological models are driven and validated for each catchment, and the best-simulation model is selected to project future streamflow scenarios. The copula function is employed to construct the joint distribution of flood peak and volume, and then the most likely realizations of bivariate quantiles are derived under different Joint Return Periods (JRPs), with the uncertainty envelope quantified with the area of 90% confidence ellipse by a copula-based parametric bootstrapping uncertainty (C-PBU) approach. Our results project an overall ascending trend of temperature and precipitation over China, and the bivariate flood quantiles and corresponding estimation uncertainty of most catchments in the future period (2056–2100) are much larger than the baseline (1961–2005), despite accompanied by substantial climate model uncertainty and spatial variability in magnitude. Many basins would be subjected to a dramatic increase of flood magnitude by over 50%, while only few basins are projected to experience a decreasing flood risk, suggesting an urgent need to increase societal resilience to a warming climate over China.