Categorization of precipitation changes in China under 1.5 °C and 3 °C global warming using the bivariate joint distribution from a multi-model perspective

This study examines the changes in the intensity and frequency of precipitation in China from a multi‐model perspective on 20 statistically downscaled fine-scale climate projections and categorizes them into four distinct patterns in response to globally targeted warming (1.5 °C and 3 °C). In a multivariate setting, the asymmetric responses of frequency and intensity to different levels of warming can be considered jointly. This study focuses on relatively moderate precipitation to determine if the ensemble of a subset of climate models, which are selected based on the categorization, can provide a better interpretation of the changing patterns compared to that from the conventional unweighted ensemble mean. The results show that the spatial distribution of the predominant category and inter-model agreement are dependent mainly on the degree of warming. As warming becomes more extensive, the projected change in precipitation tends to converge to the category that indicates an increase in both the intensity and frequency of precipitation, from the mixed-mode and even decreasing pattern. The use of subsampling to produce an ensemble of joint probability (or return period) has potential benefits in detecting asymmetric changes in the intensity and frequency of precipitation that is seen in the majority of models but hidden by the unweighted ensemble average particularly for regions where different models show mixed signals. A substantial portion of the region in China is likely to experience a transition of changes in precipitation frequency and (or) intensity under continuous warming, which would not be revealed clearly by univariate analysis.


Introduction
Reliable projections on future changes in precipitation in response to global warming are urgently needed to identify and prepare for regionally emerging challenges. Although many efforts have been made collectively to enhance the reliability of precipitation projections over the last few decades, substantial uncertainty remains in particular on the regional to local scale (IPCC 2013). From a thermodynamic perspective, the intensity of heavy precipitation will increase with increasing temperature, which is regulated by the Clausius-Clapeyron relationship that the water holding capacity of the atmosphere will increase by approximately 7% for each 1 K warming in temperature (Trenberth et al 2003, Giorgi et al 2011, Im et al 2017. While this may be true when the projected changes in precipitation are aggregated to a continental (or global scale), climate models show somewhat mixed and less coherent patterns geographically because of the overwhelming regional climate variability (Fischer et al 2014, Xie et al 2015, Lee et al 2018. More specifically, individual climate models employing different physical and dynamic processes produce a large intermodel spread in the sensitivity of precipitation to radiative forcing and the resulting change patterns. Since a particular model or response cannot be fully verifiable in the context of climate change experiments, a multi-model ensemble that equally averages all the participant models has been commonly adopted in many studies (IPCC 2013, Im et al 2017. Some studies provide evidence that the effective use of weighting may improve projections for specific applications (Giorgi and Mearns 2003, Knutti et al 2017, Eyring et al 2019. On the other hand, there is also a skeptical view due to difficulties in determining and verifying the appropriate weighting for the individual models (Weigel et al 2010). Alternatively, the subsampling ensemble has also been applied to the selection of a subset of climate projections with given certain constraints (Langenbrunner and Neelin 2017, Herger et al 2018, 2019. Most commonly, the optimized subset can be obtained by comparing the models' performance during the reference period, based on either single or multiple predictors. Despite this, it is still doubtful whether such in-sample calibration does guarantee the out-of-sample skill for the future projected simulation under climate change. Therefore, although there is an apparent consensus that characterizing the probability of future projections should be done with a multi-model ensemble approach, the methodology of how to reduce the ensemble spread remains a topic for debate. Projecting precipitation in China, which is characterized by substantial temporal and spatial variations over continental China, has been suffering from disparity among global climate models (GCMs) as well as their downscaled results. The state-of-the-art GCMs participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5, Taylor et al 2012) show limited accuracy in simulating the regional to localized precipitation and extremes in China. This is mostly because the coarse resolution of GCMs unavoidably flattens the sophisticated physiographical features of China (e.g. the sharp altitude gradient, the spatial heterogeneity of land use, and topography) (Gao et al 2006, Wu et al 2016. In this regard, dynamical or statistical downscaling methods are used to produce the fine-scale climate information that is assumed to resolve region-specific climate patterns better (Gao et al 2013, Li et al 2020. On the other hand, the individual GCMs and their downscaling results often show stark differences in the direction and magnitude of the precipitation changes in response to global warming, entailing considerable uncertainty (Gao et al 2008, Zou and Zhou 2013, Niu et al 2015. The disagreement among models appears to be greater if the precipitation intensity is rather moderate or if anthropogenic forcing (e.g. degree of warming) is not strong enough. While a significant body of studies has concluded that the intensity or frequency of extreme precipitation will increase over China, particularly under the strongly warmed climate at the end of the 21st century (e.g. Gao et al 2013, Chen 2013, Xu et al 2018, Li et al 2018, a few limited studies, targeting a near-term future or less aggressive emission scenario, have demonstrated a decrease in summer precipitation, whose regional patterns are inconsistent among models (e.g. Zou andZhou 2013, Niu et al 2015).
Changes in precipitation have been conventionally analyzed using a univariate modeling framework for a given magnitude or frequency. The univariate approach, however, may not be sufficient for representing the changes in precipitation, which can be characterized better by multiple aspects (Kao and Govindaraju 2007). In any case, because the changes in the total precipitation can be decomposed into changes in its frequency and intensity, which are mutually interconnected under the same physical processes, it is more reasonable to explore the integrated response of frequency and intensity to global warming rather than to analyze only the total precipitation or the two variables separately. For example, the future design criteria for water-related infrastructure and water resources systems should be revised if the reduction of the return period of total precipitation in the future is due mostly to the predominant increase in the intensity of precipitation rather than the increase in frequency (Giorgi et al 2019). To represent the dependence structure between random variables effectively, copulas can be used as an effective tool to estimate the bivariate joint probability (Kwon andLall 2016, Sadegh et al 2017). Several exemplary studies have been made to advocate the application of copulas for analyzing hydro-meteorological extremes over China (e.g. Zhang et al 2013, Liu et al 2014, Miao et al 2016, Yin et al 2018. This study examines the changes in the spatiotemporal statistical characteristics of precipitation, based on the joint distribution of intensity and frequency using copulas and the subsampling ensemble of the categorized precipitation changes. To represent geographically diverse features better and to reduce the substantial deficiency observed in the raw output of the CMIP5 participant GCMs, the analysis is performed using the statistically downscaled daily precipitation with a 25 km horizontal resolution covering all of China (see section 2.1). The focus of this study is to explore the relatively moderate precipitation defined as the 90th percentile precipitation of all-day year-round, instead of the rare extremes that have already been documented in many other studies. Four categories of changing patterns are classified based on the combinations of frequency and intensity change direction, which in turn serve as a standard to select the subset of models for the ensemble. The main findings are derived from a comparison of changing behaviors observed in the joint distribution of the precipitation intensity and frequency responding differently to 1.5 • C and 3 • C global warming, which are representatives of the near and far future. This study will be helpful not only in understanding the interdependency of frequency and intensity varying under the different levels of global warming but also offering new insights into the statistical analysis method of precipitation change by a combination of the joint distribution and subsampling ensemble.

Multi-model climate projections statistically downscaled
The daily precipitation simulated for the reference and future periods is obtained from the National Aeronautics and Space Administration (NASA) Earth Exchange Global Daily Downscaled Projections dataset (hereafter referred to as NEX-GDDP). The NEX-GDDP provides statistically downscaled fine-scale (25 km) climate projections throughout the longterm period  concerning most of the CMIP5 participant GCM projections (NASA 2015). Although the demand and necessity of higher resolution projections have increased, the downscaling data generated by individual research groups have limitations in wider applications in terms of the number of ensemble members and access by the public. In this regard, the NEX-GDDP could be a promising source of regional and local climate projections. Bao and Wen (2017) and Chen et al (2017) demonstrated the great potential and advantage of using NEX-GDDP against CMIP5 global projections over China from different aspects. Therefore, this dataset is used as an illustrative example for exploring the ensemble change in future precipitation. Among the 21 available models from NEX-GDDP, ACCESS1-0 is excluded since its unreasonably high precipitation for the future period compared to the other models over our target region (figure S1(available at https://stacks.iop.org/ERL/15/124043/mmedia)). The 20 selected models altogether can cover a various range of warming sensitivities to the GHGs forcing (table S1).
After the models are selected, we determine the 30 year periods corresponding to 0.48 • C, 1.5 • C, and 3 • C warming against the pre-industrial period (1861-1890) based on the global mean surface temperature of the individual GCMs obtained from the CMIP5 archive. The pre-industrial period (1861-1890) and the warming level (i.e. 0.48 • C) for the reference climate are the same one used in Sylla et al (2018) and Im et al (2019). An increase of 0.48 • C in the global mean temperature is defined as the reference climate based on the observation record of 1976-2005 compared to that of 1861-1890 (Sylla et al 2018). Warming of 1.5 • C is selected since it is regarded as the symbolic target set in the 2015 Paris Climate Agreement, while 3 • C warming is expected by a majority of GCM projections forced by the RCP8.5 scenario to be reached within the 21st century. In this regard, this study focuses on the changes in precipitation induced by 1.5 • C and 3 • C global warming, which can also represent the near and far future, respectively. Because individual models have their own phases of warming in response to GHGs forcing (i.e. RCP8.5 scenario), they will reach the target warming at different times. After determining the year surpassing 0.48 • C, 1.5 • C, and 3 • C over the pre-industrial period (1861-1890) from the 30 year running mean of the global average temperature simulated by each GCM projection, that year is assigned as the central year for the 30 year period corresponding to each warming target (table S1).
Then, the daily precipitation of NEX-GDDP over China within the corresponding 0.48 • C, 1.5 • C, and 3 • C warming periods is used for the analysis, after converting all daily precipitation less than 0.1 mm d −1 to no precipitation day. For the reference period marked as 0.48 • C global warming, the 37th highest daily precipitation of each year (including dry days with no precipitation) is selected as the equivalent of the annual 90th percentile. Here, the allday percentile is used for the threshold calculation to avoid the large sensitivity of wet-day percentiles to changes in wet-day fraction in the future (Schär et al 2016). The 30 year average of the annual 90th percentile precipitations in the reference period is used as the threshold to count the exceedance of daily precipitation at each grid point. For each warming period, the number of days in each year with daily precipitation exceeding the reference threshold is counted as the 'frequency' , and the mean intensity of these days is referred to as 'intensity' in this study. Therefore, three sets of 30 year frequency and intensity representing the different levels of global warming (i.e. 0.48 • C, 1.5 • C, and 3 • C) for each ensemble members are prepared for the analysis of joint probability distribution. Note that we exclude the very dry regions in northwest China that have more than 310 dry days per year (i.e. dry days over 85% of the all-days) during the reference period in any one of the selected models, and these grid points are marked as empty on the map in the following analyses. As this study does not intend to focus on any specific region/climate zone of China, the inclusion of the very dry region can cause problems in applying the lower threshold like 85th percentile. On the other hand, no region is excluded by the statistical significance test based on the precipitation changes because these changes often do not seem to be statistically significant, particularly in the case when the forcing level is not strong enough or in the regions where precipitation is characterized by a large variability (Christiansen 2013). However, the absence of statistically significance does not necessarily mean that there is no meaningful change. One of the main endeavors of this study is to identify meaningful changes in the characteristics of moderate precipitation obscured by the inconsistent response of models when warming is not strong enough.

Bivariate joint distribution using copula
Our focus is to categorize changes in precipitation in terms of joint probabilities to account for the dependence structure between the frequency and intensity of precipitation under climate change. Here, we utilize copulas to compute the joint return period level of the intensity and frequency of the moderate precipitation, which depends on both the reference and future precipitation over China. Copulas have been widely employed in climate and hydrology as a key component for modeling multivariate variables due to their flexibility in describing mutual dependence through a variety of copula functions. More importantly, from a theoretical point of view, the copula function approach offers another attractive feature in that different marginal distributions can be fitted for each variable, apart from the modeling of the dependence structure between variables (Kwon and Lall 2016).
Assume that the intensity and frequency, obtained for a given threshold and year, are expressed as random variables X (i.e. intensity) and Y (i.e. frequency), respectively. Given that F (x) = P [X ≤ x] and G (y) = P [Y ≤ y] are the marginal cumulative distribution function of X and Y, the joint probability can be formulated using the copula function C as follows: First, Gumbel distribution is largely identified as the marginal distributions for both intensity and frequency over China, with Akaike information criterion and the Bayesian information criterion, and their parameters (i.e. location and scale parameters) are estimated by using the maximum likelihood method. Then, this study identifies the optimum copula function with the following procedure proposed by Genest and Rivest (1993), where the copula function is constructed by Kendall's rank correlation coefficient and linear correlation coefficients. The most appropriate copula for each grid is selected by comparing the log-likelihoods obtained from the four representative copula functions (i.e. Gaussian, Clayton, Frank, and Gumbel copula) under the same marginal distribution. It should be noted that we do not limit the same copula function for the reference and the future periods. In other words, the optimal copula is selected individually for the reference and future periods. The return period is the average elapsed time or mean interarrival time between occurrences of critical events (Miao et al 2016), which are usually obtained from the inverse cumulative distribution function (i.e. quantile function) for a given return period in a univariate case. In a multivariate modeling framework, the joint return periods can be defined by estimating the probability that all the variables of interest are larger than the given thresholds. More specifically, the joint return periods of the multiple attributes, i.e. the exceedance probability of precipitation frequency and intensity exceeding a specified threshold, are simultaneously estimated as equation (2) (Salvadori and Michele 2004): where T (X ≥ x, Y ≥ y) denotes the joint return period with both X and Y exceeding the specified threshold. Figure S2 is an illustrative example that the reference 30 year time series of paired frequency and intensity (red dots in figure S2(a)) is used to estimate the full joint probability density function of P [X ≤ x, Y ≤ y] (the contour lines). Similarly, the joint probability for the future period is derived using the future time series data, and the specific return levels (i.e. 10 year/50 year return level in this study) are extracted for the comparison between the reference and future period ( figure S2(b)).

Category of precipitation change and the benefit of joint distribution
A comparison of joint distributions derived from reference (0.48 • C) and warmer (1.5 • C and 3 • C) climate conditions facilitates the categorization of the changing precipitation patterns along with the combination of increase or decrease in frequency and intensity. Figure S3  . These different categorizations have important implications for the assessment of precipitation changes in response to global warming. An increase in the total precipitation could occur under category I, but it also could be possible under Categories II or IV. Therefore, it may not be sufficient to analyze single variable in a separate manner. Multivariate modeling based categorization could identify more effectively the relative contribution of intensity and frequency that lead to the changes in total precipitation. The regional patterns and the consistency of the models dominated by each category are presented in the following section. As shown in figure S3, the changes in precipitation are characterized by multiple attributes and their relationships. The key issues are how one could systematically incorporate the dependence between attributes in understanding their changes. It is intuitively recognized that in comparison to univariate analysis the joint probability of bivariate variables could be more suitable for identifying precipitation changes simultaneously affected by changes in frequency and intensity. However, it is necessary to clearly explain how different bivariate and univariate analyses are. For this purpose, we compare the return periods derived from the different assumptions that the frequency and intensity are statistically dependent or independent. More specifically, we take a particular year and grid (i.e. 114.375 • E, 22.875 • N) with a frequency of 48 d and intensity of 43 mm d −1 during the reference period. Then, the univariate return periods for the frequency and intensity are 11 years [T (X ≥ x)] and 13 years [T (Y ≥ y)], respectively. When the bivariate distribution is taken into consideration, the joint return period [T (X ≥ x, Y ≥ y)] J is 86 years, which is less than the joint return period of 143 years that is calculated simply if intensity and frequency are assumed to be independent variables as follows.
A large discrepancy between two join return periods clearly explains that simply combining two univariate analyses together when elucidating their joint risk is inappropriate. Additionally, it implies that the intensity and frequency of precipitation should not be considered statistically independent. In order to make it more general, we compare the joint return periods across all grid points. Figure S4 clearly illustrates the effect of interdependence reflected in the joint return periods when they are calculated using bivariate distribution and univariate distribution. Figure 1 presents the predominant category of precipitation change that is classified by the 10 year return period level of joint frequency and intensity under 1.5 • C and 3 • C global warming against that of the reference period. The predominant category of each grid point is determined by the highest inter-model agreement, and the corresponding model numbers are also presented. Obviously, the spatial distribution per category and inter-model agreement show high sensitivity to the degree of warming. As warming (e.g. anthropogenic forcing) becomes more extensive, the projected precipitation change tends to converge to category I, which indicates an increase in both intensity and frequency, and to show less uncertainty with higher inter-model agreement. Under 1.5 • C global warming, the regions characterized by category IV, where the increase is only observed in intensity but not frequency, occupies a considerable area of China especially for the middle region, even though its fractions is still less than that of category I. For such region with no consistency in the intensity and frequency change, the emergence of a significant change in terms of the total precipitation might be limited because one of the quantities is likely to affect the total precipitation in the opposite direction, regardless of whether the total precipitation increases or not. On the other hand, the higher level of warming mostly eliminates the change signal of categories II, III, and IV. The domination of category I under 3 • C warming suggests that the response of models to the warming above a certain level tends to converge towards the increase in both intensity and frequency.

Results
Because inter-model agreement can infer the confidence of future projections in response to a given forcing (e.g. GHGs concentration) (Raisanen 2007, IPCC 2013, it is essential to determine if the dominant category emerges in the majority of models. As with previous studies, the inter-model agreement is higher at 3.0 • C warming than at 1.5 • C warming. In regions where only a few models support the dominant category, finding optimal ensemble method matters. At least, a multi-model ensemble with equal weighting is suboptimal because diverse categories could conceal the characteristics of a precipitation change. Figure 2 further reveals the model agreement of all the categories, respectively, which can offer a clearer idea of how to interpret the uncertainty of a model response across the different categories and how to evolve the dominant model response along with the different levels of warming. As the regions dominated by category I expand greatly from 1.5 • C to 3.0 • C warming, the inter-model agreement of category I closely resembles that presented in figure 1(d). With regard to the other categories, the extent occupied by categories II and III is reduced significantly when the warming increases from 1.5 • C to 3.0 • C. Even the areas remaining with category II or III are supported by only one or two models, which causes these two categories to totally disappear from the map of the dominant category ( figure 1(c)). It reversely implies the robust pattern in response to a certain level of warming that the physical processes employed in the majority of climate models are likely to project an increase in the intensity of precipitation, despite the change in its frequency. Interestingly, although category IV also almost disappears under the 3.0 • C global warming in figure 1(c), figures 2(d) and (h) reveal that there is no significant change in the number of models supporting it between 1.5 • C and 3 • C warming. In turn, this can be ascribed to the overwhelming increase in the inter-model agreement observed in category I caused by the transition of certain number of models from category II/III to category I, as opposed to the decrease in the number of models that produce category IV under 3 • C warming. In order to ensure that a particular threshold (i.e. 90th percentile) to define the moderate precipitation in this study would not lead to a severe distortion of the result and to compensate some arbitrary behind the selection of threshold, we repeat the same analysis using the 85th and 95th percentiles as the thresholds and the same quantities of figures 1 and 2 are also presented in figures S5-14. It can be seen that the 10 year and 50 year return period levels almost yield the same result. In terms of the different precipitation percentiles, the general patterns are also qualitatively similar. This, in turn, indicates that the analysis in this study is not too sensitive to the threshold for defining the moderate precipitation, but rather robust across a wide range of climate models. However, given a higher threshold (i.e. from 85th to 95th percentile), the strong convergence towards category I is found by reducing the extent with other categories. Put differently, both frequency and intensity are more likely to increase along the warming for the more intense part in the precipitation distribution, which is also a welldocumented characteristic of the future precipitation change (see section 4).
Note that the regions with the same number of inter-model agreement do not necessarily mean that the same models are involved. Table S2 lists the percentage occupied by the four categories concerning the individual models. Consistent with the results presented in figure 2, the most common feature noted in all models is that category I comprises the largest fraction, regardless of the degree of warming. In contrast, the majority of models (18 out of 20) produce the minimum fraction in category III under 1.5 • C warming, and such a small fraction continues to decrease (to less than 5% except for four models) with further warming. Climate models show a wellorganized pattern responding similarly to different degrees of warming, in terms of the fraction per category, unless the geographical distributions are considered. However, the models do not agree well where the same category occurs, particularly under 1.5 • C global warming, as shown in figure 2.
For in-depth analysis, two regions characterized by different behavior in terms of the category changes under different levels of warming are selected as examples for further explanation, although the two regions with limited area may not represent the big country as a whole. Along the increased warming from 1.5 • C to 3.0 • C, one region experiences a transition from category IV to I (region A in figure 1), whereas the other holds the same dominant category I (region B in figure 1). Because not all models agree with the category dominated in the selected regions, the subset of models that show the same dominant category are averaged (SENS), and it is compared with the value obtained from the unweighted ensemble of all 20 models (MENS). As addressed earlier, there are no absolute criteria for taking only a particular or all the available models for obtaining the ensemble mean. While the selection of subsampling depending on the historical performance of models is popular, the model with the best performance in the present climate will not necessarily maintain the same level of fidelity in simulating the future climate (Weigel et al 2010). In this study, more focus is placed on the model response to a given forcing for subsampling ensemble instead of confining the model performance of the present climate. The number of members in SENS for each grid is determined flexibly based on inter-model agreement. Subsampling the ensemble based on the categorization is intended to avoid the canceling out of the forced signals realized by the different categories, which is important in the study of moderate precipitation, because the disagreement among models can be large especially under relatively weak anthropogenic forcing scenario. Therefore, it can enable a clear identification of the change patterns at the given forcing level.  Figure 3 presents the joint probability of the different combinations of precipitation intensity and the frequency at 10/50 year return period level over region A. For the MENS result in the reference period, the probability that the precipitation with a 12 mm d −1 intensity occurs 52 times in a year is likely to be the same as the precipitation with a 18 mm d −1 intensity occurring 39 times, which corresponds to the 50 year return period level (0.02 exceedance probability). For this region, the predominant change category experiences a transition as the warming increases. The dominant category IV under 1.5 • C warming is converted to category I under 3.0 • C warming, which indicates a change in frequency from decreasing to increasing while the intensity remains at an increasing trend throughout the warming. Such a transition is presented in the SENS but not in the MENS, since under 1.5 • C the MENS almost remains unchanged in both frequency and intensity between the reference and warmer climates. This is possibly attributed to the cancelation of the opposite signals among the models when the forcing signal is not strong. Under subsampling, the decrease in frequency under 1.5 • C warming is stressed and yields a clearer signal. The SENS not only reveals the transition of precipitation frequency between 1.5 • C and 3 • C warming, but also shows a larger difference for the two warming periods. The frequency turns from decreasing to increasing trend, while the intensity is likely to change from a very slight increase to a much larger one, both indicating the transition toward a wetter climate with more substantial warming. The difference between MENS and SENS underscores the importance of subsampling in detecting a reverse change signal agreed by the majority of models but hidden by the unweighted ensemble average, especially for the region where the change signals among the models are quite mixed.
The quantities shown in figure 4 are the same as depicted in figure 3, with the exception that it is for the region B, where category I consistently dominates for the two warming levels and a relatively high model agreement appears (over 60% model agreement under 3 • C warming). Both MENS and SENS show that the climate in region B is likely to constantly become wetter as the warming continues, with both frequency and intensity experiencing a greater increase from 1.5 • C to 3 • C warming. SENS can yield a larger magnitude of increase under 1.5 • C warming than MENS, while such difference is smaller under 3 • C warming caused by its higher model agreement. Figures 3 and 4 show that further warming might change the decreasing signal to an increasing signal or enlarge the increasing magnitude, all of which point to an increasing precipitation. Although such increasing precipitation may also be found in the conventional unweighted ensemble average, the subsampling ensemble of the joint distribution of frequency and intensity can amplify the transition trend or magnitude of each component. This is particularly significant for regions where the models provide more diverging change signals (e.g. region A), which are also very likely to be the ones that witness more complex changes in precipitation under global warming. We acknowledge that such result is expected since we are averaging over the models showing the same signal, but this study aims to point out that averaging over a large ensemble of models without carefully distinguishing the change signals provided by individual models can sometimes lead to a neglect in the transition of change along the warming course, which may assume significance for specific regions. Chen et al (2020) also adopts the models showing the same changing directions with the multi-model median for analyses as a way to ensure the confidence levels of the results.

Discussion and conclusion
This study demonstrates the regional changes in the moderate precipitation characteristics over China in response to 1.5 • C to 3 • C global warming based on the subsampling ensemble of climate models with the same category of change patterns determined by the joint bivariate distribution of intensity and frequency. The use of joint distribution enables a more precise depiction of the interdependent relationship of the two characteristics of precipitation, and thus helps us better derive the joint risk for the future. However, if warming is not strong enough, the forced signal is very likely to be obscured by the superposed internal variability, and the consistency among the models is weakened (Fischer et al 2014). In particular, it is difficult to identify the robust response of weak or moderate precipitation changes, in contrast to heavy precipitation, which is strongly controlled thermodynamically by the increased water holding capacity at higher temperatures (Allan andSoden 2008, Giorgi et al 2011). Thus, taking an ensemble subset after applying the joint distribution is used in this study to reveal the meaningful change patterns that are often covered by the inconsistency of the models' response.
The spatial distribution of the dominant categories classified with the precipitation change patterns shows that with further warming from 1.5 • C to 3 • C, the responses of the models to forcing converge towards an increase in both the intensity and frequency and towards enhanced inter-model agreement. However, the dominance of category I under 3 • C, which may be interpreted simply as 'increased precipitation' in the conventional unweighted multimodel mean approach, may hide the transition of the precipitation frequency and intensity with continuous warming. With the subsampling ensemble that only keeps the model members in line with the majority, in the near future, precipitation frequency might decrease first compared to the present climate, and the decrease may be constrained or even converted as warming becomes more influential in the far future. Such complexity in a precipitation change is actually in accordance with the nature that the response of precipitation to global warming is a counterbalance between latent heating and radiative cooling rather than a simple linear process (Fischer et al 2014, Fischer andKnutti 2016), but may be concealed without a careful optimal ensemble approach.
The transition behavior of precipitation change emerging from the convergence to category I can be explained in the context of the changes in precipitation characteristics in response to global warming, which has already been well documented in the literature. Many previous studies have consistently found that the frequency of heavy precipitation will increase with the rising temperature, whereas the frequency of light to relatively moderate precipitation will decrease (e.g. Im et al 2008, 2017, Giorgi et al 2011, 2019, Li et al 2020. Figure S15 presents the changes in the probability of future daily precipitation along the different intensities over the two selected regions where we performed an in-depth analysis in this study. Only the changes for the intensity up to 50 mm d −1 are presented because our focus is not on the changes for the rare extremes and the increasing pattern simply continues to higher intensities. In consonance with prior studies, the projections of daily precipitation used in this study clearly show the asymmetric response in higher and weaker intensity ranges. In this regard, the behavioral change between 1.5 • C and 3 • C warmings are qualitatively similar, although the gradient of decreasing and increasing probability tends to be more pronounced in 3 • C warming than in 1.5 • C warming. For the region A where category IV prevails under 1.5 • C warming, a less gradient pattern (smaller slope) makes the probability transition (from decreasing to increasing) shifted to a more intense precipitation as compared to that under 3 • C warming. Accordingly, given higher percentile as the threshold, the extent of category IV becomes reduced (e.g. 95th percentile in figure S7). On the other hand, the transition point from decreasing to increasing probability appears in the very similar intensity for the different warming levels for region B that is consistently marked as category I under the two different warming periods, despite a less gradient for the 1.5 • C warming. Taken together, it seems that the lack of model agreement on the change category seen under the less warming does not seem to imply that the models' response to warming is uncertain. Ultimately, climate models reveal a surprisingly consistent pattern converging to category I, despite somewhat diverse patterns until the warming level is strong enough and the precipitation is intense enough. This becomes evident by comparing the results across different thresholds (85th vs 90th vs 90th percentile) and different levels of warming (1.5 • C vs 3 • C). Although the nonlinear response of the precipitation intensity under the different warming levels is not a new finding, this study reiterates the systematic transition patterns along the course of warming over China.
Overall, despite the difficulties and uncertainties existing in the precipitation projections, this study, using copulas together with subsampling ensemble, provides new insight for detecting robust changes in different characteristics of precipitation from the multi-model ensemble. Although this study provides an illustrative example for exploring the ensemble change in future precipitation, the proposed method can be applied to the analysis of other climate change impacts such as compound extremes.

Data availability statement
The data that support the findings of this study are openly available at the following URL: www.nccs.nasa.gov/services/data-collections/landbased-products/nex-gddp. Climate scenarios used were from the NEX-GDDP dataset, prepared by the Climate Analytics Group and NASA Ames Research Center using the NASA Earth Exchange, and distributed by the NASA Center for Climate Simulation (NCCS). MATLAB codes used for the analysis are available from the corresponding author upon reasonable request.