In this study, a new method is developed to generate multi-site streamflow series. First, a correlated sampling method based on spectral decomposition is employed to simulate a large number of individual years of multi-site streamflow series. Then, combined with historical design flood data, flood characteristic series considering historical floods are generated. Finally, a subset of the streamflow series corresponding to these flood characteristics is selected from the pool of multi-site simulated streamflow series. The proposed model is applied to multi-site streamflow simulations in the upper Yangtze River basin. Experimental results show that the simulated data not only ensure the fundamental statistical characteristics of the observed data, but also provide good preservation of high order self-dependent, cross-dependent at finer and coarser time scales with any marginal distributions and any correlation structures. In addition, the simulated extreme situations are more reliable for risk and vulnerability analysis of multi-reservoir joint operation systems due to the consideration of the historical design flood data.