Bivariate Frequency of Meteorological Drought in the Upper Minjiang River Based on Copula Function

: Based on the Standardized Precipitation Index (SPI) and copula function, this study analyzed the meteorological drought in the upper Minjiang River basin. The Tyson polygon method is used to divide the research area into four regions based on four meteorological stations. The monthly precipitation data of four meteorological stations from 1966 to 2016 were used for the calculation of SPI. The change trend of SPI1, SPI3 and SPI12 showed the historical dry-wet evolution phenomenon of short-term humidification and long-term aridification in the study area. The major drought events in each region are counted based on SPI3. The results show that the drought lasted the longest in Maoxian region, the occurrence of minor drought events was more frequent than the other regions. Nine distribution functions are used to fit the marginal distribution of drought duration (D), severity (S) and peak (P) estimated based on SPI3, the best marginal distribution is obtained by chi-square test. Five copula functions are used to create a bivariate joint probability distribution, the best copula function is selected through AIC, the univariate and bivariate return periods were calculated. The results of this paper will help the study area to assess the drought risk. value is 0.7177. These results will help quantify the frequency of drought events of different degrees.


Introduction
Drought is a frequent natural disaster, which affects ecology, social economy, and agriculture to a large extent. The change of drought may be faster than the average climate change with global warming [1,2]. What is more serious is that due to the expansion of the scale of industry and agriculture, social and economic development, global warming and the rapid growth of the world's population, the demand for water has risen sharply. The shortage of water resources has increased, and the global drought trend is obvious [3].
Drought is usually divided into hydrological, meteorological, agricultural, and socioeconomic drought. When the precipitation is lower than the normal level for a period of time, meteorological drought will occur [4],which may affect all other types of drought, so the evaluation of meteorological drought is important [5]. Over the past few decades, different drought indexes have been developed to assess drought conditions [6,7], including Standardized Precipitation Index (SPI) [8], Standardized Runoff Index (SRI) [9], Standardized Precipitation Evaporation Index (SPEI) [10], Standardized Hydrological Index (SHI) [11], Palmer Drought Severity Index (PDSI) [12] and so on, among which the SPI and SPEI are the most widely used [5]. According to reports, if the inter-annual temperature change in a region is not so obvious, then the results of using SPI or SPEI as research indicators will not be much different [13]. Therefore, this study chooses the SPI value as the meteorological drought assessment index.
There are many advantages of using SPI, such as simple calculations and the ability to measure drought conditions on different time scales [14]. Based on the SPI value, it is easy to extract drought characteristics, such as drought severity(S), drought duration(D) and drought peak(P) [15,16]. The analysis of drought characteristics can be univariate or multivariate. Univariate method is a traditional drought frequency analysis method [17]. However, due to the strong correlation between drought characteristics, multivariate analysis can more comprehensively characterize the drought situation. The Copula function is an excellent method for evaluating the joint probability distribution of multiple variables. Its most important advantage is that it does not need to be used on the premise that the marginal distribution of a univariate is independent [18]. At present, the copula function has been used to modeling the multivariate joint distribution of drought [19,20], flood [21,22], the joint change of precipitation and flood [23] and so on in the hydrological field.
The study area is the upper Minjiang River basin(UMR). The UMR is located in Sichuan Province, China. It is a critical water source for domestic, agricultural and industrial production in the Sichuan Basin [24]. However, the UMR has a complex geographical environment and a fragile ecological environment. Some areas have a non-zonal arid valley climate. There are large areas of arid valleys in the study area, and the foehn effect is significant [25], which makes drought become an important disaster in the area. Thus, it is a urgent need to study the drought situation in UMR.
Based on the SPI and copula function, this study analyzed the meteorological drought in UMR. The study area was partitioned into several regions based on the location of four meteorological stations using the Tyson polygon method. The monthly precipitation data of four meteorological stations from 1966 to 2016 were used to calculate the SPI values, and major drought events in various regions were counted based on SPI3 values. Drought duration, severity and peak were estimated by SPI3 value. Nine distribution functions were used to fit the marginal distributions of the three drought characteristics, and the optimal marginal distribution was obtained by chi-square test. Five common copula functions were used to create a bivariate joint probability distribution based on SPI3, and the best copula function was selected through AIC. Finally, the univariate and bivariate joint return period were calculated. The results of this study are significant to the management and distribution of water resources and the prevention of drought in UMR.

Data and Study Area
The study area in this paper is the upper Minjiang River basin(UMR). The UMR is located in Sichuan Province, China. There are many tributaries and dense river networks in the basin. It is the biggest tributary of the upper Yangtze River. The UMR is located on the southeastern edge of the Qinghai-Tibet Plateau, with high mountains and deep rivers in the area, its topography is low in the southeast and high in the northwest, which is a typical alpine valley landform [26,27]. However, due to the alternate control of the south tributary of the westerly wind, the warm Indian Ocean current, and the southeast Pacific monsoon, and under the influence of the complex and diverse geographical environment, the area has formed a unique arid valley climate feature: foehn winds in the area are strong, the atmosphere is dry all year round, and the dry and wet seasons are obvious. About 70% of annual precipitation is centralized in summer, with large annual evaporation, extreme drought in winter, and serious floods and drought disasters [25]. In addition, the UMR is located in the Longmenshan fault zone, the neotectonic movement is strong, which makes the entire mountain ecosystem fragile and changeable. In general, the UMR has a complex geographical environment and a fragile ecological environment [27]. Based on such a severe situation, the UMR was selected as the study area of this article.
The UMR basin includes all areas of Songpan, Lixian and Heishui, and parts of Wenchuan and Maoxian. There are a total of five meteorological stations. Due to the lack of precipitation data in some years in Wenchuan, this paper selects the precipitation data of other four stations as the data used in this study. Figure 1 shows that the selected four meteorological stations are evenly distributed in the UMR, which is reasonable. The four meteorological stations in UMR were used to calculate the monthly precipitation data from the daily precipitation observation data from 1966 to 2016, and the monthly precipitation data were used for the calculation of SPI.

Meteorological Drought Index Spi and Drought Characteristics
The drought index is an important variable used to assess the degree of drought and extract the drought characteristics (drought duration, drought severity, drought peak, etc.). Among them, SPI is one of the most widely used drought index s, which is recommended by the World Meteorological Organization for drought monitoring [28]. SPI was proposed by Mckee [8], its calculation is based on a multi-year monthly precipitation data series. The information of SPI response on different time scales is also different [29]. In this study, the SPIProgram downloaded from the website http://drought.unl.edu/Monitoring Tools/Downloadable SPIProgram.aspx (accessed on 15 January 2021) is only used to calculate the value of SPI on 1, 3 and 12 month time scales(SPI1,SPI3 and SPI12), the drought situation in the study area was analyzed by SPI3. Table 1 lists the SPI climate classification provided by the national standards for meteorological drought levels issued by China. According to the classification in the table, this article sets the threshold for the beginning and end of the drought time as −0.5. In addition, according to the run theory proposed by Yevjevich [30], the drought characteristics based on SPI3 is extracted. This study uses drought duration, severity, and peak to analyze drought events. The three characteristics are defined as follows: 1. Drought duration (D): the duration of SPI ≤ −0.5; 2. Drought severity (S): the absolute value of the accumulated SPI value over the duration of the drought; 3. Drought peak (P): the absolute value of the minimum SPI value during the duration of the drought. Based on the preliminary identification of drought events, in order to avoid the impact of small drought events on the analysis of statistical characteristics of drought event samples, the following treatments are made for small drought events: 1. Small drought events with drought duration of only 1 month and severity less than 1 were not included in the drought event sample; 2. When the non-drought duration between two drought events is 1 unit period and the drought severity is less than −0.2, the two adjacent drought events will be merged into one drought event.

Mann-Kendall Test
The Mann-Kendall (MK) test is often used to test the changing trends of the meteorological and hydrological time series data. Its advantage is that the tested data series don't have to follow a certain distribution [31]. The MK test null hypothesis H0 is that the change trend of the data sequence X = {X1, X2,……, Xn} is not significant. When the statistical pa-rameter|Z| ≥ 1.96, the null hypothesis is rejected within the 95% confidence interval, that is, the trend of the data series is significant. When Z is positive, it means the trend is up, otherwise, it indicates a decline in the trend [32]. This paper uses the MK trend test method to check the significance of the downward or upward trend of the SPI sequences within the 95% confidence interval. The specific calculation process of the z value is as follows [33]: The formula for calculating the variance of S is: In Equation (3), n is the number of data, k is the number of repetitions, m is the number of unique numbers (the number of groups), and tk is the number of repetitions for each repetition. When n > 10, the formula for calculating the statistical parameter Z is:

Marginal Distribution
In order to establish a binary probability distribution between drought duration, severity and peak, we must first define the univariate distribution of these characteristics. Several alternative probability distributions functions are taken into consideration in this study, namely: Weibull (wbl), Normal, Log-normal (logn), Gamma (gam), Exponential (exp), Logistic (log), Log-logistic, General Extreme Value (gev), and Generalized Pareto (gpa) distribution. In this paper, the parameters of the marginal distribution are evaluated using the maximum likelihood estimation (MLE) method. Spearman (ρ) and Kendall (τ) are used to examine the correlation between different drought characteristics.

Chi-Square Test
In order to determine the best-fitting univariate marginal distribution of each characteristics, this study uses the chi-square test to estimate the best-fitting marginal distribution. The formula for calculating the chi-square value is as follows [5,34]: Among them, n is the number of the disjoint group intervals; k is the serial number of the disjoint group intervals, Ok is the number of observations in the k-th disjoint group intervals; Ek is the expected number of observations in the k-th disjoint group intervals (according to the distribution being tested). The probability distribution function with the smallest Chi-Square value is chosen as the optimal distribution function.

Copula Function
The copula concept comes from Sklar's theorem [35]. In the Copula function, the multivariate probability distribution and the univariate marginal distribution are connected by Sklar's theorem. Then based on the joint cumulative probability distribution of the marginal distribution F1(x1),F2(x2),……,Fn(xn) (the x1, x2,..., xn are random variables), copula function can be defined [5]. Suppose that x and y are two random variables with joint distributions FX, Y(x,y) and marginal distribution functions FX(x) and FY(y), according to Sklar's theorem [36], there is a Copula function C(x,y): If FX(x) and FY(y) are consecutive, this Copula is unique. On the contrary, if FX (x), FY (y) and Copula function C (x, y) are given, the above formula defines the joint distribution function of FX (x) and FY (y) [37][38][39].
Commonly used Copula functions are generally divided into five types, including Archimedean Copula, Metaelliptical Copula, Plackette Copula, mixed Copula, and empirical Copula. Since Archimedean Copula and Metaelliptical Copula functions are easy to construct and can capture dependent structures with several characteristics, they have become very attractive functions in bivariate hydrological frequency analysis [29,39]. In this paper, three commonly used Archimedean Copula (Clayton, Frank and Gumbol-Hougaard) and two commonly used Metaelliptical Copula (Gaussian and t Student Copula) were selected, and the inference function for margin (IFM) method [40] was used to estimate the parameters of copula functions, that is, first calculate the parameter values of the marginal distribution through the MLE method, and then use the obtained marginal distribution parameters to obtain the unknown parameters in the copula functions.

Function Evaluation
The fitting efficiency of the candidate Copula function is evaluated based on the Akaike Information Criterion (AIC). The smaller the value of AIC, the higher the fitting efficiency. The calculation method of AIC is as follows [6,41]: Among them, k represents the number of fitting parameters, MSE represents the mean square error of the fitted copula function relative to the empirical copula, and XC and XE are the joint distribution functions based on the parameters and the empirical copula, respectively. MLE is the maximum likelihood of the copula function. Therefore, the copula with the smallest AIC value is the optimal copula.

Return Period
Shiau and Shen [42] proposed the return period theory of drought events. When the drought characteristic is greater than the preset value, the return period can be calculated from the expected value of the drought interval and the cumulative probability distribution corresponding to the characteristics. The calculation formula is: In the formula, E(L) is the expected value of the drought interval. FD(D), FS(S), and FP(P) are the cumulative probability distributions of drought duration, severity, and peak, respectively. TD, TS, and TP are the D, S, and P recurrence period, respectively.
According to the nature of drought, univariate analysis may cause underestimation or overestimation of drought risk [37]. Drought characteristics are related random variables, so studying the joint regression period of these characteristic quantities is more helpful to the assessment of local drought risks and the management of water resources. This article will analyze the bivariate joint probability distribution. The bivariate joint return period between drought duration, drought severity, and drought peak is divided into two situations. Here, D and S are used as examples. The combination of other characteristics is the same: (1) The return period of D ≥ d and S ≥ s is expressed by TDS; (2) The return period of D ≥ d or S ≥ s is expressed by T'DS. The calculation method is as follows [6,42]:

Temporal and Spatial Trend of Drought Situation
Based on the locations of 4 meteorological stations, the ArcGIS geographic information platform was used to generate Tyson polygons, and the study area was divided into four regions. In order to explore the changes in drought trends in various regions, this paper uses the MK trend test method to calculate the Kendall trend statistics of the SPI1, SPI3, and SPI12 at each meteorological station. The results are shown in Table 2 and Figure 2. According to the statistical distribution shown in Figure 2, the UMR can be divided into two categories, Songpan and Heishui are classified as Class I region, Maoxian and Lixian are classified as Class II region. The change trend of drought index in all time scales of Class I was significantly increasing except SPI12 in the Heishui region, which was not significantly increasing. The rising trend of SPI on the 1-month and 3-month (not cross-seasonal) timescales was more significant than the SPI on the 12-month timescales. The SPI sequence of Class II region showed a general downward trend, indicating that drought events were more likely to occur in Class II regions than before. Table 2 shows that on the time scale of 1 month and 3 months, although the drought index sequence of Maoxian region shows an upward trend, its trend rate is 0 (in fact, it is a positive number very close to 0). It can be seen that the upward trend is extremely insignificant. On a 12month (cross-season) scale, the SPI series of Maoxian and Lixian have a significant downward trend. The statistical results show the historical dry-wet evolution phenomenon of humidification in short-term and drought in long-term in the UMR.  The distribution of the dry valleys in the UMR is showed in Figure 2. It can be seen that the dry valleys are distributed in all the mainstreams of the UMR in Maoxian region, and part of the mainstreams of the UMR in Heishui and Lixian regions. The length of the dry valley in Maoxian region is the longest, followed by Lixian region. Combining the calculation results of the SPI change trend, it can be seen that, relatively speaking, the Maoxian and Lixian regions where the dry valleys are more widely distributed are more likely to become drier, that is, there is a greater risk of drought.

Meteorological Drought Assessment
SPI of different time scales reflects the different cumulative effect of drought. Figure  3 shows the SPI1, SPI3, and SPI12 sequences of the four regions. Comparing the three sequences of SPI1, SPI3, and SPI12, it can be seen that the SPI with a shorter time scale (SPI1 and SPI3) is more discrete, drought events occur more frequently, which means that the SPI with a short time scale is more capable of responding to small drought events. The long-term SPI(SPI12) treats several consecutive minor drought events as one drought event, so the long-term SPI can better reflect the long-term trend of drought, and relatively speaking, drought events last longer.   It can be found from the SPI sequences( Figure 3): in the Songpan region, from 1966 to 1972, from 1978 to 1991, and from 1996 to 2008, the SPI values were mostly negative, and the SPI values in other periods were mostly positive. This means that most of the drought events occurred in the period from 1966 to 1972, from 1978 to 1991, and from 1996 to 2008. In Heishui region, SPI was mostly negative from 1966 to 1972, from 1986 to 1987, and from 1996 to 2008, and SPI was mostly positive in other periods; in the Maoxian region, the frequency of positive and negative SPI was similar from 1966 to 1975, and the SPI was mostly negative from 1985 to 1988 and from 1991 to 2010, it can be seen that the drought lasted for a long time in the Maoxian region; in the Lixian region, SPI was mostly negative from 1966 to 1969, 1978-1980, and 1997-2012, and mostly positive in other periods. Meanwhile, Figure 3 shows that drought events occurred more frequently in Lixian during 1997-2012, and the drought was more serious.
It can be seen from the comparison of SPI sequences of different time scales, compared with SPI1, SPI3 can integrate some small drought events, which is suitable for seasonal drought and can better reflect agricultural drought scenarios. Therefore, this study is mainly based on SPI3 for drought assessment.
According to the calculated SPI value, the drought characters of D, S, and P can be extracted to evaluate the drought.  Figure 4. Figure 4 shows that the drought duration of Songpan and Heishui in the 1980s and 2000s was longer than that of other decades. Maoxian region in the 1980s and 2000s had a longer drought duration, while Lixian region in the 1960s and 2000s had a longer drought duration. The drought duration of the four regions in 2010s was relatively short. Overall, the drought duration was relatively long in the 1980s and 2000s and was the shortest in the 2010s.
All historical drought events in four regions from 1966 to 2016 were analyzed as follows. According to monthly statistics, historical drought events in Songpan mostly occurred in March; historical drought events in Heishui mostly started in January, June and December; and historical drought events in Maoxian mostly started in March June and October, historical drought events in Lixian mostly started in January, February, and April. According to seasonal statistics, drought events mostly occurred in spring and winter in the Songpan region, the proportion are 29.17% and 29.17%, respectively. In the Heishui region, the proportion of drought events that occurred in winter was 36.17%. The proportion of drought events in spring and summer was 27.66% and 29.79%, respectively, in the Maoxian region. In the Lixian region, the proportion of drought events that occurred in spring and winter was 24.49% and 32.65%, respectively. Table 3 lists some of the more serious drought events. It shows that the drought duration in Maoxian and Lixian is not only longer than that in the other two regions, but the severity of major historical drought events is also stronger than that in other regions, indicating that the drought risk in Maoxian and Lixian is relatively high. According to news reports, most of the drought disasters in Sichuan in the past 20 years occurred in the 2000s, which confirms the reliability of our above analysis. On Sichuan is facing a severe drought in spring and summer, and 5.9 million people have difficulty drinking water (http://news.sohu.com/20070404/n249186525.shtml, the accessed date is 8 July 2021). Sichuan News Net-Chengdu Business Daily reported on 27 February 2010 that since 2010, the western Sichuan Plateau has been experiencing high temperatures and low precipitation, there has been a phenomenon of droughts in autumn and winter, and the mountain snow cover was nearly 50 percent less than last year, or even at the same time for many years in February 2010.
These drought disasters have brought severe impacts on the local area in many ways: 1. Impact on humans: it has caused difficulties in drinking water for humans and animals; the dry-flow area of the Minjiang River cuts off the sources of income for residents in nearby areas who feed on and wash cars along the way; frequent "dehydration" in several sections of the Minjiang River directly affects humans when it comes to urban and rural life and industrial and agricultural production that rely on the Minjiang River for water supply; 2. Impact on agriculture: the continuous drought has caused the crops grown by local residents to turn yellow and reduce production, the supply of agricultural products is insufficient, and the price rises; 3. Impact on wild and rare animals: The construction of water conservancy projects in the UMR has changed the natural properties of the runoff and has caused a serious impact on the aquatic animals and plants of the Minjiang River. The fish species in the UMR have dropped from 40 species in the 1950s to 16 species today; 4. Impact on the environment: continuous drought has reduced the capacity of the water environment, which has aggravated the water pollution of the Minjiang River and the deterioration of the water environment. These adverse effects of drought disasters in the UMR are recorded in the above-mentioned reports.
These adverse effects of drought disasters in the UMR are recorded in the above-mentioned reports.

Marginal Distribution
To explore the joint distribution of bivariate, we must first determine the marginal distribution of univariate. Calculate the value of SPI3 in each region, using Weibull (wbl), Normal, Log-normal (logn), Gamma (gam), Exponential (exp), Logistic (log), Log-logistic, General Extreme Value (gev), and Generalized Pareto (gpa) distribution functions fit the marginal distribution of D, S, and P, respectively. The chi-square goodness of fit test was used to select the optimal marginal distribution of drought duration, severity, and peak  Table 4. Table 4 illustrates that the Exponential, Log-normal, and Log-logistic distribution were selected as the best marginal distribution of the drought duration in the four regions. Log-normal and Log-logistic were chosen as the optimal marginal distributions of the drought severity characteristics in the four regions, the best marginal distributions of drought peak were Exponential, Log-normal, and Logistic. Therefore, for the characteristic of drought duration and drought severity, it is a good choice to use Log-Logistic distribution as their marginal distribution. Log-normal distribution also has good applicability for drought peak. According to the parameter values of the best marginal distribution provided in Table 4, the value of each characteristic quantity corresponding to a specific cumulative distribution probability can be easily calculated according to needs. In order to measure the correlation between the three characteristics of drought duration, drought severity, and drought peak, the Spearman (ρ) and Kendall (τ) correlation parameters between different drought characteristics were calculated. The closer the correlation coefficient is to 1, the stronger the correlation. Table 5 shows the calculation results. The calculation results indicate that the Spearman correlation coefficients of D and S are all higher than 0.851, reaching the maximum in Heishui area (0.886), and the Kendall correlation coefficients are all higher than 0.727, and reaching the maximum value (0.757) in Songpan and Heishui regions; the Spearman correlation coefficient of S and P are all higher than 0.721, reaching the maximum in Lixian region (0.864), Kendall correlation coefficients are all higher than 0.530, and also reaching the maximum in Lixian region (0.691), which shows that there is a significant correlation between the two pairs of characteristics. Although the correlation coefficient value of D and P are smaller than the other two pairs of characteristic combinations, the correlation coefficients of Songpan, Heishui, and Lixian regions have passed the significance test of α = 0.01, and the correlation coefficients of Maoxian have passed the significance test of α = 0.05, which shows that there is a significant correlation between each characteristic. Since the positive correlation between the drought characteristics and the good fitting effect of each characteristic through different distribution functions, the copula function can be used to simulate the joint probability distribution between the drought characteristics. Table 5. Correlation coefficients among drought characteristics.

Joint Distribution of Drought Characteristics
This study used five common copula functions, Clayton, Frank, Gumbol-Hougaard, Gaussian, and t Student copulas, to set up the joint distribution of the drought characteristics based on SPI3, and the AIC method is used to evaluate the best copula function.
The AIC value in Table 6 indicates the appropriateness of t Student, Gaussian, Clayton, and Frank to establish the joint distribution of D-S. Gumbol-Hougaard is not applicable to establish the joint probability distribution of D-S at all regions. The five copula functions of Clayton, Frank, Gumbol-Hougaard, Gaussian, and t Student copulas are all suitable for describing the joint probability distribution of S-P, as well as D-P. The copula function with the smallest AIC value is selected as the optimal copula function of the bivariate joint probability distribution of each region. The best copula function and the corresponding parameters are shown in Table 7. Table 7 indicates that Gaussian and Frank copula functions are the best copula functions of D-S, as well as D-P. The best copula function of S-P is Gaussian copula function. It can be found from the optimal copula functions that for the entire UMR, the Gaussian Copula function is a good choice for simulating the joint distribution of D-S, D-P, and S-P.

Frequency Analysis of Drought Characteristics in Univariate and Bivariate
According to the best marginal distribution of each characteristic selected in Chapter 3.3, this paper gives the univariate cumulative probability distribution diagram of each characteristic through calculation. As shown in Figure 5, the cumulative probability value corresponding to the specific value of the characteristic can be read. For example, when the cumulative probability P (X ≤ x) in Maoxian area is 0.8, the corresponding drought duration(D) is 5.01 months, the drought severity (S) is 5.78, and the drought peak (P) is 1.95 .
According to the optimal copula function of the bivariate joint distribution of drought characteristics selected in Chapter 3.4, the bivariate joint probability distribution based on SPI3 (including P(D ≤ d,S ≤ s), P(D ≤ d,P≤ p) and P (S ≤ s, P ≤ p)), the joint probability distribution of the bivariate drought characteristics can be read from Figure 6. For example, when the cumulative probability of D, S and P in Maoxian area is 0.8, the bivariate joint probability P (D ≤ 5.01, S ≤ 5.78) is 0.7562, and P (D ≤ 5.01, P ≤ 1.95) is 0.6764, and P (S ≤ 5.78, P ≤ 1.95) value is 0.7177. These results will help quantify the frequency of drought events of different degrees.

Return Period Analysis
Return period analysis is an important part of drought assessment. For the determination of the drought return period, the expected value of the drought interval must first be determined. According to the analysis of drought characteristic variables, the expected drought interval E(L) for the four regions of Songpan, Heishui, Maoxian, and Lixian are 12.71, 12.97, 12.98, and 12.27 months, namely 1.0590, 1.0807, 1.0816, and 1.0226 years, which are similar to 1 year. The meteorological drought interval in southwest China is reported to be mainly affected by the superposition of monsoon and drought disturbances. The drought disturbances are mainly related to the ENSO circulation, which is closely related to the interannual planetary westerly disturbance and the interannual SST disturbance at the equator [43].
First the univariate return period (T) levels are taken to be 2, 5, 10, 20, and 50a; the corresponding values of D, S, and P are calculated, respectively, and the corresponding two-dimensional Copula function values are calculated by the optimal Copula functions of different characteristic variables. According to Equations (11) and (12), the corresponding bivariate joint return periods at a given univariate return period level are calculated.
The computed values are showed in Table 8. Table 8. Return periods of joint distribution of drought characteristics.  Table 8 illustrates that the univariate return period is between the joint return period T('and' event) and T'('or' event). The bivariate joint return period T is always bigger than T', because the calculation of the return period of the 'and' event is more restrictive than that of the 'or' event. Taking the Maoxian region as an example, the 50-year return periods of univariate of D and S are both between the TDS = 88.1931a and the T'DS = 34.8903a. In addition, under the same univariate return period level, the duration of drought in Maoxian was greater than that of the other three regions, and the severity and peak of drought in Songpan were greater than those of the other three regions, which indicates that the drought duration of Maoxian lasted longer than other three regions, and the severity and peak of the drought in the Songpan is more severe than other regions. Since the optimal marginal distribution function and optimal copula function have been obtained above, and the corresponding parameters have been calculated, in addition to the return periods corresponding to the drought events of different degrees that have been calculated in Table 8, the return period corresponding to the value of a particular drought duration, drought severity, or drought peak can also be determined according to need.

Conclusions
Drought assessment is critical to water resources planning and management. This article aims to comprehensively analyze the meteorological drought in the UMR.
In this paper, the change trends of SPI in different time scales in four regions were analyzed. The results show that the SPI sequence on a short time scale is more discrete and more able to reflect small drought events. The long-term SPI can better reflect the long-term trend of drought. The UMR showed the historical dry-wet evolution of humidification in short-term and drought in long-term. By analyzing the trend of SPI at various time scales, it is found that Maoxian and Lixian regions where the dry valleys are more widely distributed are more likely to become more arid.
Based on SPI3, the duration, severity, and peak of meteorological drought were estimated, and the drought events in each region were calculated. The results showed that the drought lasted the longest in Maoxian from 1966 to 2016, which was 183 months, the droughts in Songpan, Heishui, and Lixian lasted 166, 160, and 175 months, respectively. According to the decadal statistics of the drought duration in each region, the results show that the drought duration in the study area was relatively long in the 1980s and 2000s, and the drought duration was the shortest in 2010s. Drought events in the study area mostly started in winter and spring. Compared with the statistics of notable drought events in different regions, Maoxian not only has a longer drought duration, but also has a higher severity of historical drought events. Lixian has the highest severity of drought events in history.
According to the results of the chi-square test, this study determines the optimal marginal distribution of drought characteristics from Weibull (wbl), Normal, Log-normal (logn), Gamma (gam), Exponential (exp), Logistic (log), Log-logistic, General Extreme Value (gev), and Generalized Pareto (gpa) distribution functions. For drought duration, it is a good choice to use Log-logistic distribution as its marginal distribution. Log-normal distribution also has good applicability for drought peaks. The drought severity in different regions has different optimal marginal distributions, including Exponential, Log-normal, Logistic, and Log-logistic distributions.
Due to the dependence of the drought characteristics, this study uses Clayton, Frank, Gumbol-Hougaard, Gaussian and t Student five copula functions to fit the bivariate joint distribution to present a more realistic joint distribution result. According to the AIC value, the joint distribution of drought characteristics that is most suitable to describe each region is determined. The results show that due to differences in the correlation between drought characteristics in different regions, the applicable copula functions may also be different. For example, the optimal copula functions for D-S and D-P in different regions include Gaussian and Frank copula functions. As far as the entire study area is concerned, the Gaussian copula function is a good choice for the simulation of the joint probability distribution of the D-S, D-P and S-P.
In addition, based on the optimal marginal distribution and the optimal copula function, this paper calculates the univariate return period and the bivariate joint return period of drought characteristics to reflect the frequency of drought events of different degrees.
In general, Maoxian and Lixian have a higher risk of drought than Songpan and Heishui. According to the drought indices at different scales, almost all the SPI sequences at different scales in Songpan and Heishui showed an obvious increasing trend, while the SPI12 in Maoxian and Lixian showed an obvious trend of becoming drier. Maoxian has the longest drought duration among the four historical drought events. From the perspective of drought severity, the historical drought events in Lixian were more serious than those in the other three regions. However, this does not mean that the drought disaster in Songpan and Heishui is not serious, because except for Maoxian, the drought lasted for 183 months, the historical drought duration in the other three regions is more than 160 months, and serious drought events have occurred in all regions.
In short, the results of this paper can supply effective information for the study area to assess drought risk, so as to optimize the allocation of water resources and reduce the impact of drought on the UMR in the future.

Suggestions
Due to the frequent occurrence of drought disasters in the UMR, this article puts forward some suggestions for drought disaster management.
First of all, a good ecological environment is a strong barrier against drought disasters. Aiming at the fragile ecological environment in the UMR, new drought-resistant tree species can be cultivated, and various types of plants such as arbor, shrubs, grass, and cane can be planted to build a multi-level structure of the forest system to strengthen ecological barriers. Secondly, local residents can choose to plant crops with strong drought resistance to avoid the residents' diet from being greatly affected when drought disasters occur. Based on the concept of water conservation, relevant departments of the Sichuan government can re-allocate the limited water resource input in terms of urban and rural life, industrial and agricultural production, encourage and promote residents and factories in Sichuan to take concrete water-saving measures, and improve people's water-saving awareness.
Author Contributions: T.A., T.C. and F.Q. designed this research; T.C. and F.Q. collected the data; F.Q. analyzed the data and wrote the draft. T.A. and T.C. revised the manuscript. All authors have read and agreed to the published version of the manuscript. Acknowledgments: The authors are grateful to the editors and the anonymous reviewers for their constructive comments and suggested revisions.

Conflicts of Interest:
The authors declare no conflict of interest.