Drought Risk Assessment in Central Asia Using a Probabilistic Copula Function Approach

Abstract

The aim of this research is to adopt the Standardized Precipitation Evapotranspiration Index (SPEI) with three-month timescale (SPEI-3) to analyze drought risk in Central Asia. Based on SPEI-3, a drought event is defined through Run Theory. The multidimensional Copula function based on drought risk is then comprehensively assessed through the multivariable joint probability of drought duration, drought severity, and drought peak. Results indicate as follows: (1) the climate conditions were relatively stable from 1961–1974 and 1979–1995, while they varied from 1974 to 1979 and from 1995 to 2017, during which the study areas experienced recurrent drought. (2) The drought characteristics show noticeable spatial variability, and the severity of drought is larger in the west than in the east in Central Asia; the duration of drought contrasts with the severity of drought spatially. (3) The drought risk in the three-dimensional joint distribution is similar to the analysis using the two-dimensional distributions, and the study area has gone through the process from moderate to slight and then to severe drought risk from 1961 to 2017; the return period studied in this paper was calculated to be 80% probability in about two years.

1. Introduction

Drought has disastrous impacts on individuals’ lives and the environment. It is caused by lack of precipitation over a long period; thus, the probabilities of drought are higher in arid regions. Central Asia (CA) has attracted wide interest in its aridification, and drought is a recurrent phenomenon in CA. Under climate change, high dependence on irrigated agriculture, and increasing human disturbance, this region is susceptible to drought due to irrational water distribution.

The five CA countries, Kazakhstan (KAZ), Kyrgyzstan (KGZ), Tajikistan (TJK), Turkmenistan (TKM), and Uzbekistan (UZB), comprise the main parts of CA. This region is not only fragile in its ecological environment, but also sensitive to climate change. After the disintegration of the Soviet Union, due to the lack of unified management of water resources, drought and water resources have become key problems restricting the development and regional stability in the five CA countries [1]. The World Bank’s annual report notes that drought occurred throughout southern CA in the autumn of 2000, except for northern Kazakhstan. At the same time, rainfall was 40% lower than the average level in 2000 and 2001, and river runoff was 35–40% lower than the normal level. Although there are different studies looking at drought in CA from various perspectives, which focus on the causes, characteristics and trends of drought and try to reveal their long-term impact on agriculture and ecology in CA [2,3,4], very few studies have systematically analyzed drought risk from a multivariate probabilistic perspective. Hence, it is eager to analyze drought risk using a reasonable and comprehensive method to decrease the disadvantageous impacts of droughts in the CA.

In addition, the selection of drought index needs further consideration. The Palmer Drought Severity Index (PDSI) and Standardized Precipitation Index (SPI) are the most commonly used indices in CA for drought analysis. However, PDSI is affected by its own autoregression characteristics, and SPI is according to precipitation only and does not consider evaporation. Therefore, the choice of an appropriate drought index to analyze drought characterization is important to the research of drought risk. The newly developed Standardized Precipitation Evapotranspiration Index (SPEI) [5] provides a more accurate measure of drought that is relevant to climatic conditions, especially in arid regions. The theory of Run [6] allows one to analyze the probabilistic structure of drought durations and severities, together with the SPEI.

Drought events are characterized by many factors [7]. Independent analysis of drought factors are not able to reflect the correlations between them [8]. The univariate parametric analysis of drought events may result in over- and under-estimation of associated drought risks for water resource management [9,10]. Thus, it is difficult to apply univariate parametric analysis to assess drought risk objectively and accurately. In order to overcome this limitation, many scholars [11,12,13,14] are committed to research drought risk at higher dimensional levels. Among these approaches, the Copula function has been employed in multidimensional drought risk analysis due to its excellent characteristics [15,16,17]. Ning [18] used the two-dimensional (2D) Copula to compute the joint distribution probability of drought to analyze the regional drought risk in the northwest arid area of China. However, at present, there is less research on drought risk analysis in CA using the multidimensional Copula function. Thus, this paper uses the joint distribution probability of drought duration, drought severity, and drought peak value, which together reflect drought characteristics, to analyze drought risk in the five CA countries and to provide a meaningful attempt at a comprehensive analysis of drought risk.

The target of this paper is to employ the two-dimensional (2D) and three-dimensional (3D) Copula functions to construct multidimensional joint distributions in the five CA countries and to research the drought risk characteristics based on joint probabilities and return periods. Specifically, through SPEI and Run Theory, three drought variables can be defined. The joint distribution function will be constructed by the suitable Copula function according to their applicability in drought analysis in the region, then the multidimensional joint distribution probability of drought variables will be calculated for the five CA countries from 1961 to 2017, and the drought risk will be revealed and analyzed. An analysis of aridity changes in five CA countries will be provided in this region under climate change, which provides valuable information for the administration of water in CA.

2. Materials and Data

2.1. Study Area

This study area includes KAZ, KGZ, TJK, TKM, and UZB of Central Asia, with a collective area about 4.0 × 10⁶ km² (Figure 1) [19]. CA is home to many cross-border rivers, such as the Syr and Amu Darya rivers, which also provide daily water for most of CA’s population. After the collapse of the Soviet Union, the original integrated water management system in the region was eliminated, and the upstream and downstream countries of the river frequently disagree on water resources issues. This has a tremendous impact on the ecological environment, rational allocation of water resources, water supply, and power supply. The special geographical location, extreme weather conditions, highly concentrated population, insufficient finances, and widespread poverty also make CA extremely vulnerable to the effects of drought.

2.2. Data

There is a lack of measured data in CA in general, especially in the five CA countries. This paper uses the CRU TS 4.2.1 grid data, including monthly data of precipitation and potential evapotranspiration (PET), developed by the University of East Anglia (https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.02/). There are three reasons we use CRU data. Firstly, the CRU dataset has good quality control and homogeneity testing. It is produced by different gauge datasets [20,21]. Secondly, the resolution and time series of CRU dataset are necessary for SPEI calculations and drought studies [22]. In addition, grid data have better spatial representation and continuous availability than traditional site observation data, playing a key role in the description of drought characteristics [23]. Lastly, the CRU dataset has been widely used in CA [4,24,25] and other regions [26,27,28,29,30]. Considering that the CRU dataset produced before 1960 is not good enough [31], the calculation of SPEI needs at least 30-year datasets and the amount of Copula function input data is significant. The time series selected in this paper covers January 1961 to December 2017. The PET data of CRU are computed by a variant of the Penman-Monteith equation (http://www.fao.org/docrep/X0490E/x0490e06.htm) using the gridded data.

3. Methodology

3.1. Calculation of SPEI and Drought Classification

SPEI was developed in 2010. It not only reflects the influences of temperature fluctuations on drought events in global warming environments [32], but also detects whether droughts occur and reflects their duration over multiple time scales. The SPEI has been employed in many climatology and hydrology studies, including drought variability analysis [33], climate change [34], agriculture [35], ecological systems [36], and drought monitoring systems [37]. The two most widely used methods of PET calculation in the procedure to calculate the SPEI are those of Thornthwaite [38] and FAO-Penman-Monteith (PM) [39]. The PM equation considers the effects of temperature and other related factors, so it is more reasonable than the Thornthwaite equation, which only considers temperature [40]. In addition, the potential evapotranspiration calculated using the Thornthwaite method is temperature-based and often underestimated in arid and semi-arid areas [41]. The five regions of CA in the study area are arid and semi-arid regions. Therefore, the related R code [42], which calculates PET through the PM equation, is used in this research. (http://sac.csic.es/spei). However, specific meteorological data used by the five CA countries for the PM equation cannot be directly obtained. As the calculation of the PET data in the CRU dataset is on account of the PM method [43], this paper uses the PET and monthly precipitation data of CRU to calculate SPEI. SPEI has different time scales, which reflect different drought phenomena. SPEI-3 is usually used for water deficits caused by seasonal precipitation and temperature changes [32]. According to the data characteristics and research objectives, this paper adopts the SPEI-3 time-scale for drought risk analysis in the five CA countries. Based on the geographical characteristics of CA and relevant research, drought conditions are divided into five levels in the region according to SPEI (

Table 1. The classification of SPEI.

Level	SPEI	Category
0	−0.5 to 0.5	Near Normal
1	−0.99 to −0.5	Slight Drought
2	−1.0 to −1.49	Moderate Drought
3	−1.5 to −1.99	Severe Drought
4	<−2.0	Extreme Drought

).

3.2. The Identification of Drought Variables

Run Theory is extensively employed in the identification of drought events. In this study, drought duration (Dd), drought severity (Ds), and drought peak (Dp) were selected as the drought characteristics factors for drought risk analysis. In accordance with Run Theory (Figure 2), each factor is identified as follows:

(i)

Dd represents the duration of a drought event: the count of continuous months at which the value of SPEI is below the threshold X₀.

(ii)

Ds represents the intensity of a drought event: the absolute sum of all SPEIs during the drought period.

(iii)

Dp represents the peak value of a drought event: the minimum value of SPEI during the drought period.

(iv)

Using the above definitions, drought identification is carried out by calculating the SPEI value in the 3-month scale. The threshold is set to be −1, and drought events with drought duration less than 2 months are eliminated.

3.3. Drought Risk Probability Model

The most important step before the drought risk assessment is to carry out drought identification. In this paper, based on Run Theory, the drought event’s characteristics, including Dd, Ds, and Dp are identified. In addition, the marginal distributions of three variables are fitted. In the end, the Copula function is adopted to build the drought risk probability model.

3.4. The Theory of Copula

The Copula function was initially used by Sklar [44]. The Copula function models can construct the multidimensional joint distribution using marginal distributions and correlation framework [45]. A variety of Copula functions can be used to set up multidimensional joint distribution of drought variables. The Archimedean Copula function is one of the most commonly used Copula functions in hydrology [46] and it includes Gumbel-Hougaard (GH), Clayton, and Frank functions [47]. The three forms of Copula have been common selections for correlation models owing to their performances. The correlation coefficient of the Gumbel Copula is $2-2^{1/\mathsf{\theta }}$, and the correlation coefficient of the lower tail is zero, which gives it a strong characterization ability when describing the variation law of two random variables with upper tail correlation. The Frank Copula is a symmetry correlation function in the Archimedean Copula family. The correlation coefficients of the upper and lower tails are equal and all zero. Lower tail coefficient of the Clayton Copula is $2^{1/\mathsf{\theta }}$, and the upper tail correlation coefficient is zero [48]. Therefore, the characteristics of these three Copula functions are representative and can explain the problem from different aspects.

(1)

Parameter Estimation

The nonparametric estimation method is adopted to compute the parameters of the Copula [49] in this study. This technique is primarily related to the Copula parameter (θ) (

Table 2. Symmetric Archimedean Copula functions.

Copula	Equation
Gumbel Copula	$\mathrm{C}\left(u1,u2,u3\right)=exp\left\{-{\left[{\left(-ln{u}_1\right)}^{\theta }+{\left(-ln{u}_2\right)}^{\theta }+{\left(-ln{u}_3\right)}^{\theta }\right]}^{1/\theta }\right\}$, $\mathsf{\theta }\ge 1$
Clayton Copula	$\mathrm{C}\left(u1,u2,u3\right)={\left(u{1}^{-\theta }+u{2}^{-\theta }+u{3}^{-\theta }-3\right)}^{-1/\theta }$
Frank Copula	$\mathrm{C}\left(u1,u2,u3\right)=-\frac{1}{\theta }ln\left[1+\frac{\left(e^{-\theta u1}-1\right)\left(e^{-\theta u2}-1\right)\left(e^{-\theta u3}-1\right)}{{\left(e^{-\theta }-1\right)}^2}\right]$

). The relationship between θ and τ (Kendall correlation coefficient) is represented in the following equation. After calculating τ from the measured data, the parameters of the joint distribution can be obtained accordingly:

$$\tau =1-\frac{1}{\theta }$$

(1)

(2)

Verification and Evaluation

In this paper, in order to evaluate the fitting error quantitatively and select the appropriate Copula function, the Root Mean Square Error (RMSE) [17], Nash–Sutcliffe Efficiency (NSE) coefficient, and the Akaike information criterion (AIC) were used [50]:

$$\mathrm{AIC}:\mathrm{MSE}=\frac{1}{n-1}{\displaystyle \sum }_{i=1}^n{\left(P{e}_i-P_i\right)}^2;\mathrm{AIC}=n\log \left(\mathrm{MSE}\right)+2m$$

(2)

where m is the number of model parameters, n is the number of samples, P_i represents the Copula value of consecutive observation samples, and Pe_i represents the corresponding multivariate empirical probability. AIC is a measure of the quality of the statistical model fit. For a particular Copula function, the smaller the AIC value of the objective function value, the better the Copula function simulation effect. We define:

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left[{\tilde{y}}_i-y_i\left(\theta \right)\right]}^2}{n}}$$

(3)

$$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left[{\tilde{y}}_i-y_i\left(\theta \right)\right]}^2}{{{\displaystyle \sum }}_{i=1}^n{\left[{\tilde{y}}_i-\overline{{\tilde{y}}_i}\right]}^2}$$

(4)

where $y_i\left(\theta \right)$ is the simulated two-variable joint probability value, ${\tilde{y}}_i$ represents the empirical observation, i is the serial number of the variable, and n is the total number of variables. The range of RMSE is [0, ∞), the range of Nash–Sutcliffe Efficiency coefficient (NSE) is (−∞, 1], and when RMSE is equal to 0, NSE is equal to 1 for the perfect model.

3.5. Correlation Analysis and Establishment of Marginal Distribution Function

In this paper, in order to determine whether the drought variables are correlated and felicitous for establishing the joint distribution function, the correlation between drought variables was analyzed using Kendall rank, Spearman’s rank-order, and Pearson product-moment correlation coefficients.

The MvCAT (Multivariate Copula Analysis Toolbox) is a generalized software package which employs Markov Chain Monte Carlo simulations to estimate copula parameters [51]. It is used in this study to research the dependence structure and select the optimal marginal distribution function of the drought variable. The Exponential, Generalized Pareto, Generalized Extreme Value, Gamma, Weibull, and Gaussian distribution are used to fit the parameters of Dd, Ds, and Dp using the MVCAT tool, and the parameters are estimated by the maximum likelihood approach.

3.6. Joint Probability Distribution and Drought Risk Assessment

In combination with the definition of drought event in this paper, the drought risk probability is defined as the joint probability of Dd, Ds, and Dp. Based on the MVCAT, this paper calculates the marginal distribution for each drought variable and obtains the specific parameters of the function. The fitting effect of the function and the actual data is compared using the Quantile-Quantile (Q-Q) plot. The construction of the two-dimensional Copula function is based on the univariate marginal function by using the MVCAT calculation function. Then, comparing the RMSE, NSE and AIC values of different Copula functions, the most suitable Copula function for each set of variables is selected. The construction of the 3D Copula function is also based on the previously obtained univariate marginal distribution function [52]. However, unlike the two-dimensional case, the three-dimensional Copula function first calculates the parameter p by the inverse Kendall parameter method, and selects the most suitable Copula function with the principle of minimizing the parameter p. At the same time, according to the data characteristics, the symmetrical Copula function is suitable for constructing the three-dimensional joint distribution function. The specific formula is in Table 2.

3.7. The Return Period of Drought Event

The return period refers to how long the value of a random variable occurs over a long period. The calculation of the drought return period can provide valuable information for the rational utilize of water [53]. The single-variable and univariate return period often leads to an overestimation or underestimation of the risk rate for a given event, so this paper calculates the three-variable return period:

$$T_3=\frac{E\left(L_d\right)}{P\left\{X_i\ge x_i \vee X_j\ge x_j \vee X_k\ge x_k\right\}}=\frac{E\left(L_d\right)}{1-F\left(F\left(x_i\right), F\left(x_j\right), F\left(x_k\right)\right)}$$

(5)

where $E\left(L_d\right)$ represents the mathematical expectation of the drought interval, $F\left(x_i\right), F\left(x_j\right)\mathrm{and}F\left(x_k\right)$ represent the three univariate marginal distributions, and $F\left(F\left(x_i\right),F\left(x_j\right),F\left(x_k\right)\right)$ represents the joint probability of a three-dimensional variable group.

4. Results

4.1. Changing Trend of the SPEI

In order to reflect the characteristics of multi-scale droughts, this paper calculates the SPEI for CA on various time scales (Figure 3). Certain regularities can be found from the SPEI variation in different time scales: from 1961 to 1974 and from 1979 to 1995, the SPEI value was mostly positive, and the SPEI value was mainly negative in the other periods. This indicates that the climate in CA was more stable from 1961 to 1974 and from 1979 to 1995, while from 1974 to 1979 and from 1995 to 2017, the five CA countries experienced recurrent drought. With increasing time scale, the frequency and severity of drought events begins to decrease, but the drought duration gradually becomes prolonged. For example, the drought durations identified by SPEI-9 and SPEI-12 are longer than those of SPEI-1, SPEI-3 and SPEI-6. There is frequent alternation of positive and negative values in SPEI-1, which may lead to misestimation of the occurrence of drought events. However, the SPEI values of medium- and long-term time scales such as SPEI-9 and SPEI-12 were too rough, which ignores changes within the particular time periods. Therefore, considering the actual situation in CA, the drought risk is subsequently analyzed using SPEI-3 in this study.

Furthermore, in accordance with the classification criteria of drought (Table 1), the spatial distribution of drought occurrence frequency for different degrees of severity is counted (Figure 4). When the color is closer to blue, the frequency of drought occurring in the area is lower; when the color is closer to red, the frequency of drought occurs more frequently in the area. First, the drought frequency ranged as follows: slight > moderate > severe > extreme drought. Second, there was spatial heterogeneity in the drought characteristics of the study area. The occurrence of slight drought was higher in frequency in the eastern part of Kazakhstan, and lower in frequency in Uzbekistan and Turkmenistan (Figure 4a). The frequency of moderate droughts was evenly distributed in most parts of CA, but higher in northern Kazakhstan and lower in northern Kyrgyzstan (Figure 4b). The incidence of severe drought was generally low in CA, and the spatial distribution of frequency of extreme drought varied greatly (Figure 4c). The high values of extreme drought frequency were mainly distributed in south-central Kazakhstan, Uzbekistan, and Kyrgyzstan, while the low values mainly occurred in the north and west of Kazakhstan (Figure 4d). These results show that there are differences in drought characteristics among the countries of CA, which needs to be further studied.

4.2. The Characteristics of Drought Variables

According to the calculated SPEI-3, this paper aims to analyze the degree of drought risk; thus, the threshold for Run Theory was set at −1 and drought duration was defined to be more than 2 months, so that slight and short-duration droughts were excluded. The drought variables’ spatial distribution is displayed in Figure 5. As can be seen in Figure 6, when the color is closer to blue, the value corresponding to the variable is lower; when the color is closer to red, the value corresponding to the variable is higher. Therefore, it shows that the values of Ds varied from 3.969 to 6.625 in this area, and the distribution of Ds in the whole study area mainly ranged between 4 and 5.5. The low Ds values were distributed in the eastern area and high values were mostly in the western area, indicating that the drought severity was greater in the western areas than in the eastern areas in CA (Figure 5a). The value of Dd varied from 2.667 to 4.286 months, with the low values in the eastern area, and the middle-to-high values in the western region, which demonstrates that the duration of drought events in CA is longer in the west and shorter in the east (Figure 5b). Meanwhile, it turns out that the Ds and Dd have spatial contrasts. Specifically, the Ds in the western region is high, and the duration of drought events is also longer, which indicates that the drought situation in western CA is more serious. The value of Dp in CA was between −1.589 and 1.992; moderate drought was less frequent, but severe drought happened more frequently in CA. Extreme drought occurred in central Kazakhstan (Figure 5c).

4.3. Correlation Analysis

After the three drought variables of the five countries in Central Asia were extracted, the Kendall rank, Spearman’s rank-order, and Pearson product-moment correlation coefficients were used to measure the correlation between drought variables in each country (

Table 3. Correlation coefficients for drought events (SPEI < −1, Dd > 2).

Country	Variable	Kendall	Spearman	Pearson
KAZ	Ds_Dd	0.8446	0.9461	0.9635
	Ds_Dp	0.7073	0.8685	0.7714
	Dd_Dp	0.5799	0.7292	0.6244
KGZ	Ds_Dd	0.8476	0.952	0.9321
	Ds_Dp	0.6144	0.7971	0.7821
	Dd_Dp	0.4881	0.6343	0.6104
TJK	Ds_Dd	0.7797	0.8949	0.9177
	Ds_Dp	0.6684	0.852	0.7398
	Dd_Dp	0.4774	0.6195	0.5345
TKM	Ds_Dd	0.8214	0.9348	0.9463
	Ds_Dp	0.7515	0.9092	0.8637
	Dd_Dp	0.6009	0.7465	0.7396
UZB	Ds_Dd	0.8307	0.9428	0.9392
	Ds_Dp	0.6589	0.8411	0.8223
	Dd_Dp	0.5081	0.6691	0.663

). Coefficients closer to 1 indicate stronger correlation between the two variables. The results show that the correlation coefficients among the three groups of variables in the five countries are all close to 1, which indicates that the variables in each group show strong positive correlation among the five Central Asian countries. Thus, these three drought variables are capable of establishing the joint distribution function using the Copula function, using which the drought characteristics of Central Asia will be analyzed.

4.4. Selection of the Suitable Marginal Distribution

Selecting a suitable marginal distribution function for each variable is required to effectively construct the Copula joint distribution function. In this study, the MVCAT toolbox is used to fit the drought variables, and the optimal marginal distribution functions are selected according to their ranking. The corresponding parameter values of each drought variables in each country are also obtained (

Table 4. Variables corresponding to the parameter of the distribution function.

Country	Variable	Function	Parameter	Value
KAZ	Ds	Generalized Pareto	k	−0.0224
			sigma	2.4403
			theta	1.0840
	Dd	Generalized Extreme Value	k	3.9801
			sigma	0.0447
			mu	2.0112
	Dp	Generalized Pareto	k	−0.6429
			sigma	0.9199
			theta	0.5469
KGZ	Ds	Generalized Pareto	k	−0.5616
			sigma	5.9238
			theta	1.0529
	Dd	Generalized Pareto	k	−0.3638
			sigma	3.0777
			theta	2.0000
	Dp	Weibull	A	1.5963
			B	3.7680
TJK	Ds	Inverse Gaussian	mu	4.0423
			lambda	10.1465
	Dd	Generalized Pareto	k	−0.0343
			sigma	1.8532
			theta	2.0000
	Dp	Weibull	A	1.4730
			B	3.4742
TKM	Ds	Generalized Pareto	k	−0.3614
			sigma	4.6482
			theta	1.1508
	Dd	Generalized Pareto	k	−0.3613
			sigma	3.0205
			theta	2.0000
	Dp	Generalized Pareto	k	−1.0268
			sigma	1.5666
			theta	0.5915
UZB	Ds	Generalized Pareto	k	−0.5339
			sigma	4.9857
			theta	1.3198
	Dd	Generalized Pareto	k	−1.0902
			sigma	5.4511
			theta	2.0000
	Dp	Generalized Pareto	k	−1.1119
			sigma	1.3984
			theta	0.7025

). The Generalized Pareto function shows the highest frequency of occurrence, and most of the drought variables’ marginal distribution in five CA countries can be fitted by the Generalized Pareto function. Figure 6 shows the fitting map and Quantile-Quantile (Q-Q) plots of drought variables’ marginal distribution in Tajikistan in 1961–2017 as an example.

4.5. Selection of the Suitable Copula

The Gumbel, Frank and Clayton Copula functions are chosen to make up the pairwise 2D joint distribution of Ds, Dd, and Dp; RMSE, NSE, and AIC values for each group of joint distribution are calculated as well. The RMSE, NSE, and AIC values of the Copulas are expressed in

Table 5. The selection of 2D Copula function.

Country	Variable	Copula	RMSE	NSE	AIC
KAZ	Ds_Dd	Clayton	0.9298	0.7366	−156.2280
		Frank	0.7457	0.8306	−174.3191
		Gumbel	0.7351	0.8354	−175.4919
	Ds_Dp	Clayton	0.2717	0.9741	−257.0927
		Frank	0.2663	0.9751	−258.7414
		Gumbel	0.2743	0.9736	−256.3116
	Dd_Dp	Clayton	0.7951	0.7665	−169.0555
		Frank	0.7236	0.8066	−176.7816
		Gumbel	0.7117	0.8129	−178.1461
KGZ	Ds_Dd	Clayton	1.2267	0.5461	−137.8186
		Frank	0.5826	0.8976	−200.3612
		Gumbel	0.8258	0.7943	−171.0600
	Ds_Dp	Clayton	0.7382	0.8227	−180.4771
		Frank	0.2641	0.9773	−266.8354
		Gumbel	0.3603	0.9578	−240.7421
	Dd_Dp	Clayton	0.8359	0.7642	−170.0342
		Frank	0.5547	0.8962	−204.4800
		Gumbel	0.6748	0.8464	−188.0256
TJK	Ds_Dd	Clayton	1.1952	0.6482	−166.6986
		Frank	0.5054	0.9371	−249.3316
		Gumbel	0.7851	0.8482	−207.0471
	Ds_Dp	Clayton	0.2545	0.9811	−315.2014
		Frank	0.2214	0.9857	−328.5633
		Gumbel	0.2754	0.9778	−307.6053
	Dd_Dp	Clayton	0.8561	0.7886	−198.7347
		Frank	0.5341	0.9177	−244.0251
		Gumbel	0.6778	0.8675	−221.1546
TKM	Ds_Dd	Clayton	1.1895	0.5898	−153.6809
		Frank	0.5849	0.9008	−217.5698
		Gumbel	0.7816	0.8229	−191.4803
	Ds_Dp	Clayton	0.4255	0.9499	−246.2095
		Frank	0.2889	0.9769	−281.0619
		Gumbel	0.5783	0.9074	−218.5873
	Dd_Dp	Clayton	0.9722	0.7158	−171.8403
		Frank	0.6091	0.8884	−213.9126
		Gumbel	0.8070	0.8042	−188.6022
UZB	Ds_Dd	Clayton	1.5078	0.3384	−124.4172
		Frank	0.9986	0.7098	−159.8492
		Gumbel	1.4814	0.3613	−125.9320
	Ds_Dp	Clayton	0.3246	0.9664	−256.5018
		Frank	0.2912	0.9729	−265.8340
		Gumbel	0.4650	0.9310	−225.5920
	Dd_Dp	Clayton	1.1280	0.5810	−149.3705
		Frank	0.9410	0.7084	−164.9586
		Gumbel	1.1206	0.5865	−149.9416

. For a particular Copula function, the smaller the AIC value of the objective function value, the closer of the RMSE value is to 0, and the closer of the NSE value is to 1, the better the Copula function simulates. According to this criterion, the suitable Copula functions for each group of variables in the five CA countries are selected and highlighted in bold form (Table 5). This table shows that the Frank and Gumbel Copulas are the most suitable Copula functions of the 2D joint distribution in most cases.

The 2D joint distribution probability can only represent the relationship and probability of two drought variables. However, in order to describe and analyze the drought risk more comprehensively, it is essential to set up the multi-dimensional joint distribution of the drought variables. Therefore, the three-dimensional (3D) symmetric Copula is selected to construct the correlation structure of 3D drought variables. The parameters are calculated using the inverse Kendall parameter method and listed in

Table 6. The selection of 3D Copula function.

Parameter	Gumble	Frank	Clayton
KAZ	4.078	14.38	6.155
KGZ	3.702	12.78	5.404
TJK	3.156	10.59	4.312
TKM	4.043	14.27	6.085
UZB	3.624	12.49	5.248

. Since the smaller the parameter value, the better the fit of the model, all the drought variables in the five countries use the Gumble symmetry function to construct the 3D correlation structure.

4.6. Drought Risk Probability Assessment

The joint probability of the Copula function is used as the evaluation index of drought risk and is analyzed using the 2D and 3D joint distribution probabilities of the drought variables as follows:

(1)

2D Joint Distribution Functions

Firstly, in order to more clearly depict the probability distribution for each drought-variable group, the detailed probability distribution structure of the drought variables (Ds-Dd, Ds-Dp, Dd-Dp) in Tajikistan, which is located in the upstream CA territory, were taken as the example to be displayed in Figure 7. As shown in Figure 7, the 2D joint probability (drought risk) of the drought variables in three groups were different from each other. Specifically, the Ds-Dd joint probability was distributed centrally above 0.1, the Ds-Dp joint probability was distributed between 0.01 and 0.2, and the Dd-Dp joint probability was evenly distributed between 0.01 and 0.9.

Furthermore, the comparison results of the 2D joint distribution probability of the three groups of drought variables in the five CA countries are shown in Figure 8. It reveals that the 2D joint distribution probability of Ds-Dd, Ds−Dp, and Dd-Dp in the five countries showed similar changing trends. They decreased first and increased again, indicating that the drought risk of the five countries in CA was low in the period from 1980 to 1989 and increased significantly from 1990 to 2017. The joint distribution probability value increased significantly from the mid-1980s to the end of the 1990s, and the value was much higher after 2000 compared to that in the previous period. It indicates that the drought risk in CA became larger in recent decades. Particularly, in the period from 1980–1989 to the period 1990–1999, the drought risk in Kazakhstan increased dramatically in all three 2D joint distribution probability groups, and the probability in the joint distribution of Ds and Dd reached 0.748, which was much higher than that of other countries.

(2)

3D Joint Distribution Functions

In order to obtain comprehensive drought risk information, the 3D joint distribution probability is set up with the three drought variables. Based on Ds, Dd, and Dp, the correlation structure of drought variables in the five countries is constructed using the Gumbel symmetry Copula function (Table 2), and the drought risk (3D joint probability) value is obtained (Figure 9). Further drought risk assessment is discussed as follows.

The average trend shows that the drought risk probability of the five countries declined significantly from the mid-1970s to the mid-1980s, and then increased dramatically until around 2005. The change of drought risk is also consistent with the analyses in the 2D joint distribution: the five CA countries have proceeded from moderate drought risk to slight drought risk and then to severe drought risk during the period from 1961 to 2017. From 2000–2009 to 2010–2017, the probability of drought risk in Turkmenistan changed from the highest to the lowest in the five countries, but five countries did not change much on average. It indicates that, although the drought risk in some parts of CA has increased since the 21st Century, the overall drought risk trend is relatively stable. Particularly, the drought risk probability in Kazakhstan fluctuated greatly, with the highest value of 0.649 in the 1990s and the lowest value of 0.155 in the 1980s.

4.7. Drought Event Return Period

The drought event return period is a significant component of drought risk analysis. This paper calculates the return period of all drought events as defined by three-dimensional drought variables for all four types of drought (slight, moderate, severe, and extreme) those last more than two months. The results are displayed in Figure 10. Specifically, the return period of all types of drought events defined by three-dimensional Ds, Dd, and Dp was about 80% probability in 2 years, 15% in 2–10 years, and 5% in more than 10 years.

5. Discussion

5.1. The Application of Copula Function

Many studies have studied the characteristics of drought in CA [1,2,3,4], but no one has analyzed the risk of drought in CA from the perspective of probability. This study first assessed drought risk based on the joint distribution of multiple drought variables using MVCAT, which is an applicable and useful tool to construct the 2D Copula [51,54]. This study therefore provides a meaningful reference on MVCAT for actual drought risk analysis. Using the Copula function to compute the joint distribution probability of drought variables as a regional drought risk index has been proposed in previous research [55]. In this paper, we choose three common functions of symmetrical Archimedes functions, which are used in many studies [11,12,13,14,15,16,17], but no such attempt has been applied in CA. Therefore, this paper has provided a useful attempt for the application of a Copula function to drought risk analysis in CA for risk mitigation and water resource management. In addition, we selected other Copula functions, including asymmetric Copula functions, and compared the results to explore the different expressions for joint distribution probability of drought variables, and to establish higher dimensional Copula functions to respond to the drought risk from different levels.

5.2. The SPEI and Run Theory

In this paper, based on the CRU data, three drought variables (drought duration, drought severity, and drought peak) were defined based on SPEI-3 and Run Theory. It is clear that different time scales of the SPEI lead to different drought event definitions [34,56]. SPEI−3, used here, reflects seasonal changes and can extract sufficient drought variable series to meet the requirement of Copula data size; thus, it is recommended to use SPEI−3 for joint distribution analysis by Copulas in CA. In addition, applications of different time scales of the SPEI [34] and the integration of other drought variables (such as drought interval) for the purpose of drought risk analysis can also be further studied in future work. Many studies have used Run Theory to identify the variables of drought events [6,8,57,58], but the determination of the relevant thresholds and its duration are not the same [8,59]. Long-term and low-intensity drought events can cause serious damage [60]; however, most related studies do not focus on this issue. This paper performed a preliminary optimization of Run Theory, but a more precise definition needs to be further studied.

5.3. The Comparison between our Findings and Previous Studies

Firstly, the spatial distribution of drought frequency is consistent with the results of Guo [3] and Rumer [61] in this study. On top of that, the analytical findings of this research revealed that the drought risk in Kazakhstan from 1990 to 1999 is relatively high. On the one hand, this phenomenon may be due to meteorological and socio-political conditions [1,3,19,24]. The collapse of the Soviet Union occurred in 1991, within that period. After their disintegration, the CA countries underwent great political and economic changes, and the water resources once managed and dispatched begun to compete in various countries; in addition, the problem of cross-border water broke out strongly during that time [4]. On the other hand, this phenomenon may have been influenced by the country’s boundary [62] since Kazakhstan is the largest of the five CA countries. In addition, in the selection of the marginal fitting function of the drought variable, our findings were similar to the fitting function used in the arid region of Northwest China [63]. Since the arid region of Northwest China and the study area of this article belong to central Asia, the climate and environmental conditions are very close. Therefore, this reflects that the results of this article are more reasonable. Finally, in the selection of Copula functions, the three functions we have selected have also been used in arid regions [14]. The drought risk in the three-dimensional joint distribution is similar to the analysis using the two-dimensional distributions, which also verified the rationality of the selection of the 3D Copula function in this study.

6. Conclusions

In this paper, based on CRU data, three drought variables (Ds, Dd and Dp) were defined based on SPEI-3 and Run Theory. Their joint probability was constructed using the Copula function as the evaluation index of drought risk in the five CA countries. The drought risk in CA in the past half century is analyzed, and the findings are as follows:

(1)

By calculating the 1, 3, 6, 9, and 12-month scale of SPEI, it is found that climatic conditions were relatively stable during 1961–1974 and 1979–1995, while they varied more from 1974 to 1979 and from 1995 to 2017, during which the five CA countries experienced recurrent drought. With the increase of the time scale, the frequency and severity of drought events began to decrease, but the drought duration gradually became prolonged.

(2)

The severity of drought in CA is greater in the west than in the east, and the duration of drought is spatially contrasted with the severity of drought. In areas with high (low) drought severity, the drought duration is also high (low). The drought events in CA are mainly severe; moderate droughts are less common, and extreme droughts mainly occurred in central Kazakhstan.

(3)

For drought risk analysis based on multi-dimensional joint probability of drought variables, the drought risk in the three-dimensional joint distribution is similar with the analyses using the two-dimensional joint distribution. The five CA countries have gone from moderate drought risk to slight drought risk and then to severe drought risk from 1961 to 2017. Furthermore, the return period of drought events, defined by three-dimensional Ds, Dd, and Dp, was calculated at about 80% probability in 2 years, 15% of 2–10 years, and 5% for more than 10 years.