Abstract Machine learning (ML) models have been widely used for hydrological simulation. Previous studies have reported that conventional ML models fail to accurately simulate extreme flows which are crucial for design flood estimation and associated risk analysis in the context of climate change. Therefore, this study proposes a joint probabilistic rainfall‐runoff model (JPRR) for improving high‐to‐extreme flow projection. With the aid of paired copula constructions, bootstrap aggregation, and multi‐model ensemble approaches, the proposed model is able to effectively characterize the dependence relationships of predictors (i.e., precipitation time series with different moving sums) with various probability distributions. JPRR has been applied to four pristine basins in China, representing different climate zones and landscapes. The results reveal that JPRR significantly outperforms three well‐known ML models (i.e., random forest, artificial neural networks, and long short‐term memory) in high‐to‐extreme flow simulations. In JPRR, the copulas exhibiting the right tail dependence play a more important role in streamflow simulations at mountainous basins. Moreover, a significant difference in streamflow projections (from 2030 to 2099) derived from JPRR and benchmark models imply that flood risks from conventional ML models may be underestimated under changing climatic conditions.

Abstract What elements should a parsimonious model reproduce at a single scale to precisely simulate rainfall at many scales? We posit these elements are: (a) the probability of dry and linear correlation structure of the wet/dry sequence as a proxy reproducing the distribution of wet/dry spells, and (b) the marginal distribution of nonzero rainfall and its correlation structure. We build a two‐state rainfall model, the CoSMoS‐2s, that explicitly reproduces these elements and is easily applicable at any timescale. Additionally, the paper: (a) introduces the Generalized Exponential ( ) distribution system comprising six flexible distributions with desired properties to describe nonzero rainfall and facilitate time series generation; (b) extends the CoSMoS framework to allow simulations with negative correlations; (c) simplifies the generation of binary sequences with any correlation structure by analytical approximations; (d) introduces the rank‐based CoSMoS‐2s that preserves Spearman's correlations, has an analytical formulation, and is also applicable for infinite variance time series, (e) introduces the copula‐based CoSMoS‐2s enabling intermittent times series generation with nonzero values having the dependence structure of any desired copula, and (f) offers conceptual generalizations for rainfall modeling and beyond, with specific ideas for future improvements and extensions. The CoSMoS‐2s is tested using four long hourly rainfall records; the simulations reproduce rainfall properties at multiple scales including the wet/dry spells, probability of dry, characteristics of nonzero rainfall, and the behavior of extremes.