Modeling Joint Relationship and Design Scenarios Between Precipitation, Surface Temperature, and Atmospheric Precipitable Water Over Mainland China

This study employs copula model to discuss joint relationship and design scenarios between precipitation (P), surface temperature (T), and precipitable water (PW) in seven regions and mainland China over 1980–2016 within a bivariate framework. Markov Chain Monte Carlo modeling in a Bayesian framework was applied to calculate copula parameters while best fitted copulas were assessed using Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Root Mean Squared Error (RMSE), Nash‐Sutcliffe Efficiency (NSE), and maximum likelihood criteria. Results showed that the spatial variations of P decreased from South to Northwest, T decreased from South to North and from East to West. Distributions of P and T were similar, with higher values in regions c and d. PW was higher in the south, with more than 64.54 mm in the southwest, 25 mm in the north‐central. The correlation between PW and P as well as PW and T were higher than 0.74. Birnbaum Saunders distribution was considered appropriate to fit P in regions a and f while P and PW in other regions are good with Generalized Pareto. T performed better with generalized extreme in all regions. Moreover, Gumbel, Frank, and Joe copulas provided good fit for PW and P in regions b, f, and g, respectively while in regions a, c, d, and e, Nelson provided perfect fit. For PW and T, Nelson was a good fit in all regions. The bivariate probabilities and design scenarios in various regions suggest tremendous variations, with regions having high probabilities associated with low return periods. We showed that joint analysis of climate variables gives more robust design quantiles.

when all the vapor in the overlying air column is assumed to condense under favorable condition (Trenberth, 1998). The analysis of PW has been achieved in many ways including satellite data inversion, water vapor mutation, and atmospheric cloud physics (Gong et al., 2017;Li et al., 2010;Ma, 2007), which are all aimed at establishing a numerical relationship between water vapor content and P as well as the relation between PW and water vapor pressure (Fan et al., 2016;Wang & Liu, 1997;Yang et al., 2010) in typical climate processes.
Also, Cai et al. (2004) showed that the PW across Tibet and its neighboring regions differed significantly with seasons while Wang et al. (2006) and Gong et al. (2017) investigated the seasonal pattern and progression of PW in the Qilian mountain area and believed that atmospheric water resources in spring and summer have great development potential. Therefore, PW could be valuable to characterize the degree of climate dryness and wetness (Fan et al., 2016;Liu et al., 2013;Yang et al., 2008), and also reveal processes of atmospheric water vapor cycle and atmospheric water resources utilization (Fan et al., 2016;Gong et al., 2017;Zhang et al., 2014). Consequently, P, T, and PW directly affect water resources availability and perform a dominant function in any climatic units as a major response factor with moisture dynamics and radiative impacts (Trenberth et al., 2005).
Despite this progress, there is a need for a reliable and accurate characterization of the joint relationship and dependencies between these variables through an in-depth analysis within the theoretical frameworks (Leonard et al., 2014). Notably, univariate methods may not be adequate to describe the design return periods of interdependent variates (Grimaldi & Serinaldi, 2006;Salvadori et al., 2016), therefore, multivariate methods are often used to understand their associations. The copula concept, as one of the major approaches of multivariate analysis, has gained popularity and have been proven valuable in describing and analyzing the structure of dependences between hydrology and climatology variables (De Michele et al., 2005;De Michele & Salvadori, 2003;Favre et al., 2004;Salvadori & De Michele, 2004a, 2004b. Although several copula families can be used to express various probabilistic characteristics of variables, however, most studies have employed few, which includes Clayton, Frank, Gaussian, Gumbel, t-copulas among others. Copula classes vary in their complexities and dependence structure, having parameters ranging from one to three. According to Balakrishna and Lai (2009), parameters of copula indicate the strength of the joint relationship among variables and the most commonly used parameter estimation method includes: local optimization, such as the Newton-Raphson method (Gräler et al., 2013;Salvadori & De Michele, 2004b), L-moment-based approach (Brahimi et al., 2015), global optimization and Bayesian analysis (Kwon & Lall, 2016;Parent et al., 2014), which performs better than the traditional methods with regards to bias and computing costs.
Copulas models have been broadly employed in several applications, which includes rainfall modeling and evaluation (Li et al., 2013;Vernieuwe et al., 2015), multi-index drought evaluation (Hao & AghaKouchak, 2013), multivariate frequency analysis (Parent et al., 2014), extreme value analysis, and coastal flooding risk evaluation (Jongman et al., 2014). Therefore, since copula has been extensively employed to model the joint connection between two or many variables, our primary intent in this study is to accurately characterize the joint relationship and dependencies between P, T, and PW climatic variables using copula models. In some cases, some variables may possess dependence structures that are asymmetrically skewed, which cannot be represented by the commonly used copulas; therefore, we employed a broad range of copula families having differing degrees of freedom that could take the dependences of symmetric and asymmetric variables into consideration.
Specifically, we (i) select appropriate marginal distributions for P, T, and PW, (ii) review 20 copulas families including those that have not been widely used in climatological fields, (iii) estimate parameters of copulas utilizing Markov Chain Monte Carlo (MCMC) modeling in a Bayesian framework, (iv) rank the copula models in order of performance using AIC, BIC, RMSE, and NSE goodness of fit measures, (v) use the best fitted copula to discuss the joint relationship and design scenarios between characteristics of P, T, and PW in various climatic regions and entire mainland China over 1980-2016 within a bivariate framework. The remaining parts of this paper are organized thus: section 2 briefly explains the study area and data sets, marginal distribution selection, copula concept, and multivariate dependence analysis. We AYANTOBO ET AL. 10.1029/2020EA001513 2 of 26 defined the Bayesian analysis and MCMC simulation algorithm for copula parameter estimation. Next, we highlighted the goodness of fit measures, and finally introduce a methodology for the joint probability and return periods. We continue with the study results in section 3 which discusses the joint relationship and dependencies between P, T, and PW in terms of probabilities and return periods. Section 4 ends this article with a summary of the joint relationship modeling between the climatic variables and presents some conclusions.

Study Area and Data Sets Processing
The climatic regions of China fluctuate significantly due to its huge space and intricate topography (Zhai et al., 2010a). Western and Northern climates are dry while Eastern China is semihumid and humid (Wu et al., 2011). According to climate and topography, mainland China is separated into seven regions which includes Northeast humid/semihumid warm (region a), North China humid/semihumid temperate (region b), Central and Southern China humid subtropical (region c), South China humid tropical (region d), Steppe zone of Inner Mongolia (region e), Northwest desert (region f), and Qinghai-Tibet Plateau (region g) (Ayantobo et al., 2017(Ayantobo et al., , 2019Zhao, 1983). Regions a-d and e-f are called the Eastern Monsoon and Northwestern regions, individually ( Figure 1). P and T fluctuations are obvious since the greater part of China's regions are situated in the East Asian monsoon climatic belt (Wu et al., 2011).
The Climate Forecast System Reanalysis (CFSR) product was given by NCEP, which consolidated significant improvements like a higher resolution of T382 (∼38 km) having 64 vertical heights spanning from surface to 0.26 hPа. Based on a grid point interpolation scheme, CFSR uses NCEP coupled forecast model, comprising а modular ocean and spесtrаl aerial model through 3D variational data assimilation (Saha et al., 2010(Saha et al., , 2014   Considering the duration of data availability, 1980-2016 periods are selected and the variates selected for evaluation are daily precipitation, air temperature, and specific humidity. In the actual atmosphere, because the water vapor content in the atmosphere above 300 hPa is little, so 300 hPa is used as the upper boundary of the atmosphere. Specific humidity data of eight isobaric surfaces (1,000,925,850,700,600,500,400,300 hPa) are superposed to get the PW. Since reanalysis provides data sets four times daily relating to 00:00, 06:00, 12:00, and 18:00 UTC, the daily mean P (mm/day), air T (K/day), and specific humidity (kg kg −1 ) are calculated by averaging the subdaily values at these four times daily. More information about the reanalysis used in this study is shown in Table 1.

Estimation of Precipitable Water (PW)
The PW is described as the total water vapor in the air column on the ground, which represents the depth of water quantity that could be obtained when all the water vapor in the air column is assumed to condense or fall. It is an essential index to characterize atmospheric water resources and the degree of climatic dryness and wetness. PW can be expressed as (Fan et al., 2016;Liu et al., 2013;Yang et al., 2008;Zhai & Eskridge, 1997): where PW is the water vapor content of the whole atmosphere per unit area (mm); q is the specific humidity (kg kg −1 ); Pt is the atmospheric top pressure (i.e., 300), Ps is 1,000 hPa, and g is the acceleration of gravity (m/s 2 ). According to the calculation reference method in Zhuo et al. (2013) and Fan et al. (2016), this study calculates monthly PW.

Description of Marginal Distribution Function
The probability and design return period estimation generally requires fitting marginal distribution functions. As against many studies that use few marginal distributions (Zheng et al., 2015), employing a broad range of distributions is necessary to reduce initial assumptions on data through opting for the best fitted distribution function. Here, we use 15 varieties of continuous marginal distribution functions to determine the most appropriate one that maximally fits P, T, and PW variables. The summary mathematical description of the distribution functions with their corresponding parameter range are given in Table 2.
Readers are referred to Johnson et al. (1993Johnson et al. ( , 1994 for a review and mathematical formulations of these distributions. The marginal distributions parameters are calculated using the algorithm of maximum likelihood which reduces the gap between empirical probabilities and the modeled equivalents. Utilizing the Gringorten position-plotting method (Gringorten, 1963), the empirical nonexceedance probabilities are compared against the fitted distributions as follows:  and  are shape, scale, and location parameter, where k ≠ 0 , , and k are shape, scale, and location parameter. f or < x, when k > 0, or for  < x <  - /k when k < 0 9 Exponential  Here, P k is the probabilities of kth value, k is kth smallest value in the data ordered sequentially, and n is the size of data sets. The best marginal distributions are chosen using the BIC metric, while a visual contrast of the empirical probability versus fitted distribution, and QQ plots, were adopted to confirm how acceptable are the distribution fits.
Likewise, a joint cumulative distribution, H, can be used to construct copula as C(u 1 , u 2 ) = H[F 1 −1 (u 1 ), F 2 −1 (u 2 )]given [u 1 , u 2 ] = [F 1 (x 1 ), F 2 (x 2 )]. The description of copula can also be extended to d-variables as: In this study, we employed 20 bivariate copulas to build models for P, T, and PW. A summary description of the bivariate copulas illustrating various forms of dependence structures are given in Table 3 (Joe, 2014;Nelsen, 2007). Here, we maximally adapt the computed joint probabilities to their corresponding empirical values. Hence, this approach is a nonparametric way of fitting copulas that employ pseudo-observations to derive parameters of copula.

MCMC Modeling in a Bayesian Framework
The Bayesian investigation has widely been applied in various disciplines (Ellison, 2004;, including hydrology Kavetski et al., 2006;Thiemann et al., 2001), for model inference goals. This theorem connects all modeling uncertainties to the parameters and calculates posterior distribution of the model parameters via: AYANTOBO ET AL.
10.1029/2020EA001513 6 of 26 x a b x f x a b e a a a and b are scale and shape parameters, where a > 0 and b > 0  Li et al. (2013) 13 Gumbel (1960) and Barnett (1980) 15 is coined evidence. Evidence, a fixed value can be removed during analysis if the aim is to calculate parameter posterior distribution, and posterior parameter distributions are calculated as: A flat uniform prior could be employed in case there is a lack of data on the prior distribution of parameters (Thiemann et al., 2001). If there is no correlation in error residuals or if the Gaussian has zero mean, the likelihood function could be evaluated as (Sorooshian & Dracup, 1980): Here  is the standard deviation error estimation. For numerical stability, equation (6) can be logarithmically translated into: where  can be calculated using equation (8) considering underlying assumptions regarding error residuals.
Equation (7) can then be further simplified as: The constants can be removed and the log-likelihood function given as: Because of the difficulties in solving Bayes' equation (5) analytically, the Markov Chain Monte Carlo (MCMC) simulation is adopted to sample from the posterior distribution, hence solving the Bayes' equation numerically. The MCMC is a kind of statistical approach used to inspect from high dimensional composite distributions (Andrieu & Thoms, 2008;Sadegh et al., 2017). The simulation optimally finds an estimation of the global optimal at little computational account. For a detailed description of MCMC simulation algorithms refer to Sadegh et al. (2017Sadegh et al. ( , 2018 Huynh et al. (2014) 17 Nelsen (2007) 18 Plackett (1965) 19

Evaluating Copula Models Performances
In this study, we employed AIC, BIC, RMSE, NSE, and maximum likelihood value to examine the performances of the various copula models. The likelihood value is estimated using equation (6), where the maximum likelihood reduces the residuals between observation data and model simulations. AIC, in contrary to the likelihood value, presents a further robust measure of the quality of model forecasts, and is expressed as (Aho et al., 2014;Akaike, 1974Akaike, , 1998: Here D is the parameter number of the model while ℓ is value of log-likelihood of the best parameter set in equation (7). By employing the Gaussian theory of error residuals in equation (8) and a constant CS, the equation can be simplified as: Similarly, according to Schwarz (1978), BIC is formulated as: and can be rewritten into NSE and RMSE have also been widely used and they focus on the minimization of residuals as follows: These metrics appraises copula performances regarding how close empirical observed probabilities (Y) are to their modeled bivariate counterparts (  Y ). Metrics with the smallest AIC, BIC, RMSE, and a close to 1 NSE value is considered a more suitable model fit. Immediately after selecting the best copula for each region, then the joint occurrence relationships are then built.

Bivariate Joint Occurrence Probability
In this study, two kinds of occurrence bivariate joint probabilities (Shiau, 2006) for PW and P pair are evaluated and described as: In the same vein, the joint probabilities between PW and T are estimated according to equations (17)-(20)

Univariate Frequency Analysis of Events
The return period is an analytical model of the anticipated recurrence interval of a particular event across an extensive period. This statistical theory has been generally employed for infrastructure design and risk analysis plans. Accordingly, the univariate return period is expressed as (Ayantobo et al., 2017;Sadegh et al., 2018;Shiau, 2006): In which T PW , T P , and T T represents the univariate return period and are described as a function of the E(L d )and marginal distribution probability of PW, P, and T, respectively. E[L d ] is defined as the expected interarrival period of observed cases (Salvadori et al., 2011(Salvadori et al., , 2014. The concept of a univariate design return period could be expanded to a much higher dimension frequency investigation by using copula functions to consider the dependence between variables.

Multivariate Frequency Analysis of Events
In bivariate frequency analysis, each variables' marginal distribution is evaluated empirically, and the best fitted copula is calibrated on the empirical joint probabilities to create the joint density functions, thereby modeling the dependence structure. The bivariate return periods are obtained and defined for two cases: T PW&P and T PWorP , and are expressed as (AghaKouchak, 2015;Salas et al., 2005;Shiau, 2006):   Figure 2 shows the distinct spatial patterns of mean monthly P, T, and PW based on CFSR reanalysis over mainland China during 1980-2016. Generally, P distribution decreased from the South toward the Northwest, and the highest was observed over the Southwestern, which was also the maximum P regions in China (Dai et al., 2006). The maximum mean monthly P was 17.44 mm while the least P was recorded in the Northwestern parts, where the mean monthly P was less than 1.00 mm.

Spatial Distribution of Mean Monthly P, T, and PW
The complicated topography of China influenced T distribution and it was shown that region g having an altitude of about 8,806 m had lowest T (i.e., around 259.52 K) than region d (i.e., about 298.07 K), located at about −154 m. This showed that T decreases from South to North and from East to West of China. The distributions of P and T were similar and the values were regularly higher in Southern China (regions c and d) than the other regions.
As latitude increases, PW decreases with a wide variation between Southern and Northern China, with more than 64.54 mm in the Southwest, 25 mm in the North-central region. The Southwest PW has shown an overall increasing trend, while the decrease of PW over Northern China has been related to a decreasing P trend over the past decade (Ayantobo et al., 2020). Due to the influence of the East Asian monsoon, studies have also shown that P and PW changes over China exhibited distinct seasonal variations, with summer PW significantly higher than other seasons (Ayantobo et al., 2020). Apart from the impact of T in the distribution of PW over eastern China, altitude also plays a significant role. The spatial distribution of PW was consistent over regions a, b, e, and f, where the altitudes were low with minimum PW values of less than 12 mm. This confirms that altitude and topography affect PW distribution, which has an important influence on the thickness and moisture content of the air column (Liang et al., 2006). P, T, and PW are important climate parameters, and the study of changes in their spatiotemporal pattern is crucial due to climate change (Cong & Brady, 2012;Pandey et al., 2018) and the knowledge about their interrelationship is vital for planning natural ecosystems. Knowing the variation of T and its relationship with other climate variables is very important for climate planning (Balyani et al., 2017). The P, T, and PW have very high temporal variations (Balyani et al., 2017), and the dependence of P on T causes its high variation over the time (Cong & Brady, 2012;Pandey et al., 2018). The dependence between these climatic parameters makes their analysis problematic especially when they are treated as independent of each other (Cong & Brady, 2012;Pandey et al., 2018).
In this study, our attention is focused on the joint relationship between P, T, and PW in the climatic regions of China within the copula framework in order to appraise their dependencies which is critical to reliable modeling and prediction. By employing a multivariate framework approach, variables of P, T, and PW are combined using their joint distribution (Hao & AghaKouchak, 2013). Therefore, we followed this structure to build the joint occurrence probability distribution of P, T, and PW, which can assist in the risk assessment of extreme hydrological and meteorological events. In the next sections, we explained how we empirically evaluate the marginal distribution of individual variable and then build joint distributions of P and PW, and T and PW using copula functions.

Dependence Analysis Between P, T, and PW
Modeling correlation structure between P, T, and PW are vital for the proper characterization of atmospheric water resource dynamics. Here, we investigated the interdependency of these three variables by various AYANTOBO ET AL.

Table 4 Values of Kendall, Spearman and Pearson's Correlation Coefficients for Precipitable Water (PW; mm) Versus Precipitation (P; mm), and Precipitable Water (PW; mm) Versus Temperature (T; Kelvin) Variables Over Regions a-g and Entire Mainland China (h) Between 1980 and 2016
correlation techniques namely, Spearman's rank, Pearson, and Kendall's rank correlation coefficient, each of which displayed meaningful interdependence among the variables. As shown in Table 4, the Kendall, Spearman, and Pearson correlation coefficient demonstrated the links amid the various variable combinations. Generally, the values recorded by spearman were the highest while Kendall was the lowest. However, the correlations between pairs of PW and P, and PW and T combinations are significantly higher than 0.74 in all the regions including the entire mainland China.
The PW and T were highly connected with coefficients varying from 0.959 (region g) to 0.984 (region a). These values are relatively higher than the correlation between PW and P, which ranged from 0.915 (region a) to 0.945 (region g). This agrees with previous studies where a significant correlation between water vapor content and P was observed (Fan et al., 2016;Zhai & Eskridge, 1997). It should be recorded that P variations may change PW via surface evaporation, an outcome that might be critical in the dry areas of Northern China. The PW and P changes may be the outcome of the same atmospheric dynamics (Fan et al., 2016;Zhai & Eskridge, 1997). Overall, the results revealed that an increase or a decrease in P is connected with an increase or a decrease in PW. Similarly, an increase or a decrease in T was connected with a corresponding increase or decrease in PW in most regions.
AYANTOBO ET AL.  Table 5 The

Marginal Distribution Selection for P, T, and PW
Estimating the occurrence probability or design return period of an episode generally entails fitting distribution functions. So, we fitted 15 kinds of continuous marginal distributions functions to P, T, and PW and selected the best fit using the BIC metric. The marginal distributions parameters are evaluated via a AYANTOBO ET AL.

10.1029/2020EA001513
13 of 26  maximum likelihood procedure that reduces the gap separating empirical probabilities and their modeled equivalents. Table 5 shows the best fit marginal distributions for P, T, and PW with their calculated parameters over regions a-g and entire mainland China from 1980 to 2016. It could be observed from our study that Birnbaum Saunders distributions were considered good for fitting P in regions a and f while P and PW in other regions including entire mainland China are best fitted with Generalized Pareto.
Concerning T, Generalized Extreme values are relatively uniform and performed well for all regions including entire mainland China. The QQ visual plotting in Figures 3a-3c as well as the plots of the best fitted distribution versus empirical probability amounts in Figures S1-S3 (SI), were adopted to check the acceptableness of the distribution fit for P, T, PW, respectively in regions a-g. It could be seen that the theoretical (red line) and empirical probabilities (blue dots) were pretty close and their relations were near the 45° diagonal line. These marginal distributions were later used for bivariate frequency analysis between the variables.

Selection of Copula Models
The bivariate copulas and regional marginal distributions in Tables 2 and 4, sequentially, were employed to ascertain the best fit copula for creating occurrence bivariate joint distribution between P, T, and PW anomalies. We evaluated 20 bivariate copula models and the copula parameters were estimated using MCMC simulation in a Bayesian framework. One copula was chosen for each region using the performance metrics procedure discussed in section 2.6. Details of the selected copulas including the copula parameter, p value as well as associated RMSE, NSE, MAX_L, AIC, and BIC metrics for PW and P, and PW and T in each region and entire mainland China are provided in Table 6. It was worth noting that after inspection of 20 copula families fitted separately to PW and P, and PW and T, only four copulas (Nelson, Gumbel, Frank, and Joe) are considered the best model for describing the asymmetric dependence structure, through the consistency of all the goodness of fit criteria.
Indeed, the Gumbel (RMSE = 0.12, NSE = 1.00), Frank (RMSE = 0.14, NSE = 1.00), and Joe (RMSE = 0.16, NSE = 1.00) provided a very good fit for modeling PW and P in region b, f, and g, respectively, and were thus selected as the best copula. In regions a, c, d, and e, it was noticeable that Nelson provided a perfect fit for modeling PW and P. The Nelson copula parameter yielded an RMSE and NSE value of 0.12 and 1.00, 0.13 AYANTOBO ET AL.

Table 6 Best Fitted Copula Families for Precipitable Water (PW; mm) Versus Precipitation (P; mm) and Precipitable Water (PW; mm) Versus Temperature (T; Kelvin) Over Regions a-g and Entire Mainland China (h) Between 1980 and 2016
and 1.00, 0.25 and 1.00, 0.11 and 1.00, respectively in these regions. For PW and T, Nelson provided a good fit in all climatic regions, where the NSE values were all equal to 1.00 (NSE = 1 is connected with an ideal fit) and the associated RMSE ranged from 0.09 (region b) to 0.18 (region g). It was worthy of note that in entire mainland China, Nelson copula seems to be a robust and adequate for simulating PW and P, and PW and T variables. The Nelson copula had a p value, RMSE, and NSE of 0.01, 0.15,1.00 and 0.05, 0.13, 1.00 for PW and P, and PW and T, respectively. An interesting observation from our study was that different copula families with similar performance metrics can yield various dependence structures in modeling PW, P, and T.
Based on statistical and graphical comparability of 19 bivariate Archimedean copula functions, Ayantobo et al. (2018) also found that Nelson copulas provided a good fit for drought duration and severity in virtually all the climatic regions of mainland China. On the contrary, Aghakouchak et al. (2010) employed t and Gaussian copula models of Elliptical class for the bivariate analysis of T and P, and discovered that t copula was far better than Gaussian. Similarly, Cong and Brady (2012) illustrated five copula models including Gumbel, t, Normal, Clayton, and Frank as best models for determining the dependence among T and P variates. Using these copula functions, Pandey et al. (2018) modeled the structure of dependence between T (Minimum, Maximum, Mean) and monthly P, of which the Normal distribution was the best fit, so was employed for bivariate modeling of T and P variates. All these studies emphasize the significance of the choice of copula in any dependence analysis. Hence, in the rest of this study, the dependence structures of PW, P, and T are then modeled employing the chosen copula in individual regions to get associated occurrence joint probabilities and design return periods, which is imperative for assessing and forecasting regional scenarios.

10.1029/2020EA001513
16 of 26 . Bivariate joint probability isolines between precipitable water (PW; mm) and precipitation (P; mm) (i.e. P PW&P (PW ≥ pw, P ≥ p)) derived from bestfitted copula over regions a-g and entire Mainland China (h) between 1980 and 2016. The joint probability isolines are color coded with joint density levels with blue representing lower densities and red denoting higher densities. Black dots are observed pairs of PW and P.

Figure 4b
. Bivariate joint probability isolines between precipitable water (PW; mm) and precipitation (P; mm) (i.e., PWorp ( ≥ or ≥ )) derived from best fitted copulas over regions a-g and entire Mainland China (h) between 1980 and 2016. The joint probability isolines are color coded with joint density levels with blue representing lower densities and red denoting higher densities. Black dots are observed pairs of PW and P.

10.1029/2020EA001513
18 of 26 Figure 5a. Bivariate joint probability isolines between precipitable water (PW; mm) and temperature (T; Kelvin) (i.e., PW&T ( ≥ , ≥ )) derived from best fitted copula over regions a-g and entire Mainland China (h) between 1980 and 2016. The joint probability isolines are color coded with joint density levels with blue representing lower densities and red denoting higher densities. Black dots are observed pairs of PW and T. Figure 5b. Bivariate joint probability isolines between precipitable water (PW; mm) and temperature (T; Kelvin) (i.e., P PWorT (PW ≥ pw or T ≥ t)) derived from best-fitted copula over regions a-g and entire Mainland China (h) between 1980 and 2016. The joint probability isolines are color coded with joint density levels with blue representing lower densities and red denoting higher densities. Blue dots are observed pairs of PW and T.

Regional Joint Probability of Events
To estimate the bivariate joint probabilities in different regions, we applied equations (18) and (20) which was based on the best fitted copula family selected for modeling the dependence structures between PW and P as well as PW and T discussed in section 3.4. The bivariate occurrence probabilities are represented with contour lines because various PW, P, and T composite could result in different probabilities. The bivariate joint probability isolines derived from the best fitted copulas associated with PW and P for   and 5b for regions a-g between 1980 and 2016. Color code is used to describe the joint probability isolines, with joint density levels of blue and red signifying lower and higher densities, respectively. A scatter plot is included in the plot to show the level of correlation between the variables.
It could be seen that the joint occurrence probabilities varied with regions, with P or significantly higher than P and values. Applying the probability isolines plots outcomes for a particular scenario, there are several relations of PW and P as well as PW and T. For example, considering region g in Figures 4a and 4b where we used the Joe copula to model PW and P design scenario, and analyze the corresponding probabilities. Based on these figures, two joint probabilities were associated when the pw and p exceed 30 and 6 mm respectively. The P PW&P was 0.03 ( Figure 4a) while P PWorP was 0.50 ( Figure 4b). On the other hand, when the Nelson copula model was used to model PW and T in the same region, the probabilities associated with pw and t exceeding 30 mm and 285 K, respectively for P PW&T was 0.20 ( Figure 5a) and that for P PWorT was 0.48 (Figure 5b). The significance of these bivariate occurrence joint probabilities has also been recorded by many studies (Ayantobo et al., 2018;Liu & Zhang, 2011;Shiau, 2006) because the joint probabilities are indeed valuable in giving valuable guidance for evaluating water resources systems. Awareness of these various joint probabilities are very important for the management and assessment of risks imposed by extreme meteorological and hydrological events, management of productivity of agricultural products, establishment of an early warning system for extreme events such as flood and drought, and also useful for given valuable guidance for evaluating as well as improving water resources systems condition in future.

10.1029/2020EA001513
19 of 26 T PWorP case joint return period isolines derived from Nelson copula for entire mainland China. Both joint return period cases are color coded with joint density levels with blue representing lower densities and red denoting higher densities. Black dots are observed pairs of PW and P.

Regional Univariate and Bivariate Return Period of Events
Here, based on the best fitted copula family selected for illustrating the dependence structures between PW and P as well as PW and T, we propose a bivariate expected event with respect to different return period levels (see equations (22) and (23)) in different regions of China. The likely scenarios quantiles refer to events that are not extreme and that is usually expected to happen in any considered year, and are importantly less than the main extreme events. The bivariate return period isolines for 5, 10, 20, 50, and 100 years derived from the best fitted copulas associated with PW and P for T PW pw P p Utilizing the contour plots outcomes for a particular return period, there are different associations of PW and P as well as PW and T. Figures 6 and 7 shows the univariate and their associated bivariate joint return period isolines of PW and P, and PW and T pairs, respectively, derived from Nelson copula for entire mainland of China based on the distinct procedures detailed in section 2.7. The color codes describe the design return period isolines with blue and red signifying lower and higher density levels, respectively. A scatter plot is included in the plot to show the level of correlation between the climate variables. According to 6.40 mm, and T; 294.39 K) highlighting that ignoring the interactions between PW and P and also PW and T could lead to underestimation and overestimation of the events. These differences for a 25 years design return period although not very significant but cannot be overlooked. This shows that copula function allows providing the dependence structures between two variables.
AYANTOBO ET AL.    can we explain the high PW over regions c, d and most especially g. PW is the basis of P formation and it could be observed that their distribution changes could be influenced by many factors, such as moisture flux, monsoon, temperature, surface evaporation (Zhai & Eskridge, 1997;Zhou et al., 2009), atmospheric dynamic conditions for P, unstable energy (Cai et al., 2004;Fan et al., 2016;Zhai & Eskridge, 1997), and elevation (Cai et al., 2004).

Implications for Univariate and Bivariate Design Scenarios
For instance, region g is the Qinghai-Tibet Plateau region, and Wu and Zhang (1998) observed that the Plateau continuously draws wet air from the low latitude Oceans in summer. During winter, а tenacious anticyclone which occurs over the Arabian Sea, near the shore of Somalia, carries water vapor through the westerly stream toward the Plateau region. Besides, the water stored within the Plateau itself could enter the water cycle and goes into the atmosphere (Ayantobo et al., 2020;Wu & Zhang, 1998). Moreover, since the mean global surface T has risen by 0.2-0.3°C across the last 40 years with a significant warming trend over the Plateau since the 1980s (Zhang et al., 2013), we may infer that the elevated water vapor across the Plateau could be because of increment in surface air T across the Plateau, which enables the air over the Plateau to retain more extra moisture. The implication of these is that the climate becomes humid over the Plateau, which makes PW more abundant and increases P potentials. Therefore, P, T, and PW play a dominant role in climatic systems as a major response factor with radiative impacts and moisture dynamics (Trenberth et al., 2005), as it directly affects water resources availability.

Summary and Conclusion
Climatic variables are interdependent and knowing their joint dependencies is usually requisite to a more reliable modeling and forecast. Since copula has been used broadly for building the joint relation between two or more variables, in this paper, we used 20 copulas families including those that have not been widely used in climatological fields to discuss the joint relationship and design scenarios between characteristics of P, T, and PW in mainland China over 1980-2016 within a bivariate framework. Our study employed the MCMC algorithms in a Bayesian framework to evaluate copula parameters, and the copula performances in each climatic region were assessed using maximum likelihood, AIC, BIC, RMSE, and NSE criteria. We later used the best fitted bivariate copula to establish a joint relationship and analyze frequencies of P, T, and PW in different climatic regions of mainland China. Here we summarize our conclusions.
The spatial variations of P, T, and PW were different from region to region. Generally, P decreased from the South toward the Northwest, and the highest can be observed over the Southwestern. The complicated topography of China influenced T distribution, making T decreasing from Southern to Northern China as well as from East to West. The distributions of P and T were similar and the values were generally higher in regions c and d. The PW contents are higher in the South than in the North, with more than 64.54 mm in the Southwest, 25 mm in the North-central region.
The correlation coefficients between pairs of PW and P as well as PW and T are significantly higher than 0.74 in all the regions, with PW and T values (0.959 (region g) to 0.984 (region a)) relatively higher than PW and P (0.915 (region a) to 0.945 (region g)). The implication is that an increase (decrease) in P is connected with an increase (decrease) of PW, while an increase (decrease) in T are connected with an increase (decrease) of PW in most regions.
The selection of marginal distribution performs an influential role in defining design quantiles. Birnbaum Saunders distribution is considered appropriate to fit P in regions a and f while P and PW variables in all other regions including the entire mainland China are best fitted with Generalized Pareto. Concerning T, the generalized extreme values performed well for all the regions including entire mainland China.
Copula choice is very critical in multivariate event assessment. It was worth noting that after examination of the 20 copula families fitted separately to PW and P, and PW and T, only four copula families (Nelson, AYANTOBO ET AL.

10.1029/2020EA001513
22 of 26 Gumbel, Frank, and Joe) are considered the best model. Gumbel, Frank, and Joe copula provided a very good fit for modeling PW and P in regions b, f, and g, respectively while in regions a, c, d, and e, Nelson provided a perfect fit. In the case of PW and T, Nelson provided a good fit in all the climatic regions of China. We suggest employing a broad range of models with various features to guarantee the fitted multivariate model is representative.
The bivariate occurrence probabilities and design return periods in various regions suggest huge disparities as regions having high probabilities are oftentimes connected with short return periods and vice versa.
One key question from our study is that to what extent the selection of marginal distribution and copula model affects the estimated probabilities and design return periods. We noted that few copula classes with somewhat comparable performance metrics may give significantly diverse kinds of probability and return periods isolines from one climatic region to another. This proposes the issue of which model should be entrusted to independently present the design values. This study has shown that these discrepancies could originate from various causes including errors from the goodness of fit of both marginal distribution and copula model, model structural errors, and also the observed joint probability errors (Sadegh et al., 2018). Currently, most studies examine quite a few copula models in their investigation, which may result in huge discrepancies in predicted joint probabilities and design return periods values. Such discrepancies could be reduced by rigorous marginal and copula fitting, taken from a broad range of copula families, to ensure that best fitted ones are selected for modeling. We show that a copula analysis gives a more advanced design quantile and therefore should be embraced by researchers. Ultimately, these estimates could be transparently reported to those liable for assessment to provide a reference role in the design evaluation of national and regional water resources as well as developing strategies in order to reduce the risk of climatic and hydrological phenomena in future.