Meteorological drought analysis using copula theory and drought indicators under climate change scenarios (RCP)

The study was carried out to assess meteorological drought on the basis of the standardized precipitation index (SPI) and standardized precipitation evapotranspiration index (SPEI) evaluated in future climate scenarios. Yazd province, located in an arid region in the centre of Iran, was chosen for analysis. The study area has just one synoptic station with a long‐term record (56 years). The impact of climate change on future drought was examined by using the CanESM2 of the CMIP5 model under three scenarios, that is, representative concentration pathways RCP2.6, RCP4.5 and RCP8.5. Given that a drought is defined by several dependent variables, the evaluation of this phenomenon should be based on a multivariate analysis. For this purpose, two main characteristics of drought (severity and duration) were extracted by run theory in a past (1961–2016) and future (2017–2100) period based on the SPI and SPEI, and studied using copula theory. Three functions, that is, Frank, Gaussian and Gumbel copula, were selected to fit with drought severity and duration. The results of the bivariate analysis using copula showed that, according to both indicators, the study area will experience droughts with greater severity and duration in future as compared with the historical period, and the drought represented by the SPEI is more severe than that associated with the SPI. Also, drought simulated using the RCP8.5 scenario was more severe than when using the other two scenarios. Finally, droughts with a longer return period will become more frequent in future.


| INTRODUCTION
On a global scale, droughts have increased in recent decades, but this trend is more pronounced in arid and semi-arid regions (Dai, 2011;Saadat et al., 2013). Drought is an extreme event and is a creeping phenomenon as compared with other natural disasters, which greatly affects the environment and human life. A drought is characterized by severity, duration and frequency. These characteristics are not independent of each other, so that drought is a complex natural disaster (Mishra and Singh, 2010). The effects of drought, as compared with other natural disasters, will have a wider range (Wilhite et al., 2014). Therefore, in drought analysis, it is necessary to consider its multivariate nature and spatial variability explicitly.
The prediction of drought conditions in future is considered a vital step in preventing future damage. The global warming of recent years, which has been consolidated by several researchers through climate change model analysis (Solomon et al., 2007), will lead to an increase in the occurrence of droughts (Sheffield and Wood, 2008). The determination of drought characteristics is very difficult because several drought indicators are available for evaluation, so that introducing a unique index is very difficult (Heim, 2002). The standardized precipitation index (SPI) and standardized precipitation evapotranspiration index (SPEI) are the most commonly used indicators for analysing the severity, duration and frequency of meteorological drought. The SPI is more accessible than other indices (Caccamo et al., 2011) because it only uses precipitation in the definition of drought. Also, at different time scales it identifies different drought types (Vicente-Serrano et al., 2010). The short time scales are related to soil moisture, which indicates agricultural drought, the medium time scales indicate hydrological drought and the long-time scales indicate groundwater drought (Edwards and McKee, 1997). The SPEI, in addition to precipitation, uses the temperature parameter to identify drought periods; it also includes the effects of temperature variability on drought assessment (Dubrovsky et al., 2008). Therefore, it can consider the effect of climate change on drought better than the SPI (Vicente-Serrano et al., 2010).
The future climate situation at the global scale is simulated by general circulation models (GCM) based on the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (AR5), which was completed in 2014 (IPCC, 2014). According to the AR5 models, the increase of surface temperature at the end of the 21st century is expected to be > 1.5 C for most scenarios. In the preparation of the AR5, which was gradually released between 2013 and 2014, the output of the CMIP5 models was used. These models use new emission scenarios called representative concentration pathways (RCPs) (Moss et al., 2010). These scenarios have four pathways: RCP2.6, RCP4.5, RCP6 and RCP8.5, which are named according to their radiative forcing in 2100 (Van Vuuren et al., 2011). The GCM is the most appropriate method for investigating the effects of climate change on environmental systems at global scales (Sehgal et al., 2018). However, the application of these models is not suitable at regional or local scales owing to their coarse grid spacing (Hay et al., 2000;Gebremeskel et al., 2005), which is now of the order of 10 km (Cavicchia et al., 2014).
The most important and suitable tool for linking the local/regional scale and GCM scale is downscaling. Several methods have been proposed for downscaling that consider the temporal and spatial mismatch between the global and regional/local scales. In general, two statistical and dynamic methods are considered for this purpose. In the statistical method, the empirical/statistical relationship between large-scale and local/regional variables is established (Wetterhall et al., 2006;Palatella et al., 2010). Among statistical downscaling methods, the SDSM has been widely used for the downscaling of climate variables all over the world (Huang et al., 2011).
As mentioned above, drought characteristics are dependent on each other; therefore, drought analysis based on one of these characteristics as a univariate cannot describe the dependence structure between them. Therefore, the simultaneous assessment of drought characteristics due to its multivariate nature is necessary and the results of the multivariate analysis of drought characteristics are important for assessing its potential risks (Shiau and Modarres, 2009). Given that the drought characteristics may have different marginal distribution functions, the most suitable tool for achieving the joint behaviour of the drought characteristics is multivariate copula functions. Copula functions as a multivariate analysis method are widely used to analyse the functional structure of drought characteristics (Mirakbari et al., 2010;Chen et al., 2013;Dodangeh et al., 2017;Hao et al., 2017).
In the present study, meteorological droughts are analysed by the SPI and SPEI in the present and future periods in Yazd province of Iran. In order to assess the drought characteristics in the present, precipitation and temperature data of Yazd synoptic station were extracted in the period 1961-2016. The CanESM2 (Second Generation Canadian Earth System Model) data from the AR5 were downscaled to examine drought changes in future (2017-2100). Bivariate copula functions were used for joint modelling of meteorological drought based on the SPI and SPEI indices. By the joint probabilities derived from the copula functions, the bivariate and conditional drought return period were calculated based on its characteristics in order to investigate the future drought changes in comparison with the historical period.

| Study area
Yazd province, with the area of 129,285 km 2 , is situated on an oasis where the Dasht-e Kavirdesert and the Dasht-e Lut deserts meet ( Figure 1). The mean annual precipitation is 57 mm for the 63-year period (1954-2016); also, the mean annual temperature is 19.58 C; the maximum temperature in the warmest month is 45.6 C; the minimum temperature in the coldest month is 16.0 C; and the maximum wind speed is up to 120 kmÁ hr −1 . The average annual relative humidity is 25-30%. The mean annual evaporation is 2,500-3,500 mm, which is higher than precipitation. The higher evaporation rate compared with precipitation has caused the study area to be classified as an arid region.
Owing to the geographical and climatic conditions, the study area has been affected by prolonged and severe droughts in recent years (Dastorani et al., 2011). According to Negaresh and Khosravi (2011), the mean precipitation in the study area has decreased by 7.5% over the past 50 years and mean temperature has increased by 1.0 C in this period. The drought studies in Yazd province have used the SPI for drought analysis as a univariate (Dastorani et al., 2011;Shayegh and Soltani, 2011;Nassimi and Mohammadi, 2016;Asadi, 2017;Javanmard et al., 2017). However, owing to the multivariate nature of drought the univariate methods for assessment of its characteristics do not fully describe the condition of drought. Based on the literature review, little research has been done to study the effects of climate change on drought in Yazd province. Among the studies, only Dastorani et al. (2011) examined the impact of climate change on meteorological drought by two SPI and RDI indicators based on the HadCM3 model derived from the 4th Assessment Report in Yazd province (Solomon et al., 2007). In the present study, precipitation, minimum and maximum temperature were projected by HadCM3 model under A2 and B2 scenarios for the period 2010-2039. The results of the present study showed the drought has an increasing trend under the A2 scenario, and a decreasing trend under the B2 scenario compared with the reference period . In general, no single study was found that uses copula theory for meteorological drought analysis in Yazd province under new climate change scenarios.

| Data
To predict the trend of drought under climate change conditions, three types of inputs were used. These include daily precipitation and temperature observed data in the period 1961-2016, National Centre for Environmental Prediction (NCEP) reanalysis data in the period 1961-2005, and CanESM2 model data from the AR5 in the period 1961-2005 for the historical and 2006-2100 for the future period. Observational data (daily precipitation and temperature, then converted to monthly scale) were extracted from Yazd synoptic station for 56 years and the simulated data were obtained by downscaling of the CanESM2 model under three RCP scenarios using the SDSM. Drought periods were extracted from the SPI and SPEI in the two periods. Finally, future drought changes in comparison with historical ones were investigated based on the two characteristics, severity and duration, using copula theory. Figure 2 describes the flowchart of the methodology.

| Downscaling of the CanESM2 model data using the SDSM under the RCP scenarios
In order to produce precipitation and temperature data for future drought analysis, the CanESM2 model data were downscaled using the SDSM. (The CanESM2 is the Second Generation Canadian Earth System Model provided by the Canadian Centre for Climate Modeling & Analysis-CCCma.) This model has three scenarios: RCP2.6, RCP4.5 and RCP8.5. The SDSM is used to create quantitative relationship between large-scale variables/ GCM and small-scale variables (local/regional scale) (Wilby et al., 2002). The method has four main parts, which include determination of the NCEP predictor variables, calibration and verification of model, and finally simulation of precipitation and temperature data under the RCP scenarios for the future period.

| Calculating the SPI
The SPI is used to identify and evaluate drought events at multiple time scales (McKee et al., 1993). It is calculated on the basis of the monthly precipitation deficit relative to its long-term average (Seiler et al., 2002). To calculate this index, the Gamma distribution function was fitted to the monthly precipitation time series, and the cumulative probability of precipitation was converted into the standard normal distribution with mean zero and standard deviation (SD) of 1 (Equations 1 and 2). The gamma distribution function is selected as the best fit to the precipitation data (McKee et al., 1993;Angelidis et al., 2012;Stagge et al., 2015;Yacoub and Tayfur, 2017): where H(x) = q + (1 − q) G(x), where q is the probability of zero precipitation; and G(x) is the cumulative probability derived from the gamma distribution. Also, the constants are C 0 = 2.515517, C 1 = 0.802853, C 2 = 0.010328, d 1 = 1.432788, d 2 = 0.18929 and d 3 = 0.001308.

| Calculating the SEPI
Based on the SPI, precipitation is the only parameter that influences meteorological drought. This is based on the assumption that the variation in precipitation is higher than other climatic parameters (Vicente-Serrano et al., 2012). However, climate change studies showed that an increase in temperature has a significant effect on the global water cycle and the frequency of dry periods. The SPEI (Vicente-Serrano et al., 2010), in addition to precipitation, uses an average temperature parameter to determine drought conditions. It is calculated at different time scales based on the non-exceedance probability of the water balance (D), which means D results from the difference between precipitation (P) and potential evapotranspiration (PET) (D = P -PET) (Equation 4): where w = (−ln (P)) 0.5 p ≤ 0.5; and p is the probability of exceeding a determined D value. The constant values are same as the SPI. The SPI and SPEI values are categorized into dry, wet and normal periods according to McKee et al. (1993). The most important advantages of the two SPI and SPEI methodologies are the flexibility in choosing time scales. In the present study, three time scales, six, nine and 12 months, were used to extract drought periods from these two indicators.

| Definition of drought characteristics based on run theory
A drought event based on the SPI and SPEI is considered to be a period in which the values of the indices are continuously negative and ≤ −1 (Nedealcov et al., 2015;Lee et al., 2017) that is described by severity, duration and frequency. Run theory is a probabilistic method (Yevjevich, 1967) used to extract drought characteristics based on the SPI and SPEI. The method is widely used to calculate the severity, duration and frequency of drought while determining a threshold level ( Figure 3). According to run theory, drought severity is equal to the sum of the SPI/SPEI values below the threshold level; the length of time during which the SPI/SPEI value is continuously below the threshold level is drought duration; and the frequency of drought F I G U R E 2 Flowchart of the methodology is the number of events when the SPI/SPEI is below the threshold level.

| Bivariate drought analysis by copula theory
For multivariate probabilistic modelling, Sklar (1959) proposed copula theory. The copula function makes it possible to combine several univariate distributions from different families to create a bi-or multivariate distribution, taking into account the dependence between variables. In other words, the copula function (C(u 1 , u 2 , …, u N )) is a connection function for the relationship between the distribution functions of the random variables X 1 , X 2 , …, X N with the marginal distribution functions Fx 1 (x 1 ), Fx 2 (x 2 ), …, F XN (x N ), defined by Equation (5) (Nelsen, 2006): Many types of multivariate distributions model the dependence structure of variables on the assumption that the marginal function is uniform. This assumption causes errors in multivariate analyses. The most important advantage of using copula functions is that the dependence structure of variables is made in the absence of the same kind of marginal distribution functions of the variables, indicating that the creation of a joint probability distribution function does not require the uniform marginal functions of each variable. By copula theory, marginal distribution functions are chosen to create multivariate distributions, as well as to describe nonlinear and asymmetric relationships between variables. Copula functions have several families, including the Elliptical (t-copula, Normal), Archimedean (Gumbel, Clayton, Frank, Ali-Mikhail-Haq), Extreme Value (Husler-Reiss, Galambos, Tawn, and T-EV, Gumbel), and other families (Plackett, Farlie-Gumbel-Morgenstern) (Yan, 2007). Of these, the Archimedean and Elliptical families are most commonly used (Madadgar and Moradkhani, 2013). The Archimedean copulas have symmetric and asymmetric forms. The symmetric forms refer to one-parameter copula, and asymmetric forms have more than two parameters. The Elliptical copulas do not have closed form in contrast to Archimedean copulas. The advantage of Elliptical families is that this group includes correlation in both the upper and lower tails. In the present study, the symmetric Archimedean copulas (Frank, Gumbel, Clyton, Joe) and Elliptical copula (Gaussian) were used for bivariate modelling of meteorological drought. When a copula function is used to create multivariate distribution, a correlation between variables is required. In the absence of correlation, it is not possible to use the copula function. In the present study, Kendal's correlation co-efficient was used to determine the correlation between drought severity and duration.

| Estimation of the copula parameter
To estimate the copula parameter, two non-parametric and parametric methods were used. In the non-parametric, the relationship between the generator function of each copula (ϕ) and the Kendall correlation co-efficient (τ) (Equation 6) were used (Genest and Rivest, 1993). In the parametric method, the parameter was estimated using the maximum log-likelihood function (Equation 7) (Favre et al., 2004). In these equations, c θ is the copula density function; F is the marginal distribution function; and X 1k , X 2k , …, X pk (k = 1, …, n) are the random variables: 2.3.7 | Selection of the copula function To choose the most suitable copula function, which has the best fit to the drought severity and duration, the F I G U R E 3 Definition of the drought characteristics based on run theory empirical joint probability of the two variables was calculated by the empirical copula function (Equation 8), and was compared with the theoretical values obtained from Archimedean and Elliptical copulas. In order to compare the empirical and theoretical copulas, the normalized root mean square error (NRMSE) (Equation 9), Akaike information criterion (AIC) and Bayesian information criterion (BIC) (Equations 10 and 11) (Akaike, 1974;Bozdogan, 2000) were calculated. In these equations, S and D represent the empirical probability of drought severity and duration (respectively, P ei is the empirical copula; P i is theoretical copula); k is the number of model parameters; n is the number of observations; and L is the maximum log-likelihood function: where 1 is the indicator function that takes value equal to 1 when its argument condition is satisfied.
2.3.8 | Bivariate analysis of the return period of drought A drought event may occur several times throughout the year and continue for several months. The return period is the average interval between occurrence of a phenomenon, which is repeated several times over time and used to describe the severity and frequency of drought (Liu et al., 2015). The bivariate analysis of return period of drought is based on the average interval time between the drought events (E(L)) and the joint cumulative distribution function of the drought severity and duration (C (S, D). Based on the definition of the return period, if N is the length of an event, n is the number of events and L is the interval time between events, so E(L) is the average interval time between successive drought events (E(L) = N/n). Thus, bivariate analysis of drought return period is calculated according to Equations (12) and (13) (Shiau, 2006). When the two variables, severity and duration of drought, are considered simultaneously, the bivariate return period is based on both the severity and duration that exceed a certain value (S ≥ s and D ≥ d) (Equation 12), or the return period, where the severity or duration exceeds a certain value (S ≥ s or D ≥ d) (Equation 13): 3 | RESULTS

| Downscaling of precipitation and average temperature data
Considering that modelling the climate change requires a long-term record and it is also required at least 30 years' of monthly precipitation to calculate the SPI (Svoboda et al., 2012), the authors just used Yazd synoptic station over the study area , which has an acceptably long record. After reviewing and controlling the quality of the observation data, the NCEP predictor variables, which had the highest correlation with each observation data, were selected. The number of correlated NCEP variables depends on the length of record and the type of process (conditional and unconditional) ( Table 1). The calibration and verification of the model were based on the NCEP variables. In this step, in order to increase the accuracy of the model, different variance inflation factors for precipitation and temperature parameters, and bias correction for precipitation only, were selected. The variance inflation factor controls the magnitude of the variance inflation in downscaled data. It changes the variance to regression model estimates of the local process (SDSM manual). The SDSM implements bias correction and variance inflation techniques to reduce the standard error (SE) and to increase the amount of variance in order to obtain the best result of downscaling (Dibike and Coulibaly, 2005). Based on the results of the validation of the model, the best value for each factor was determined by statistical comparisons. Thus, the variance inflation factors for temperature and precipitation were 6 and 10, respectively, and the bias correction for precipitation data was considered as 1. Figure 4 shows the performance of the SDSM for precipitation and temperature data during the historical period . The SDSM performance for temperature is higher than precipitation because of the continuity of the temperature variable and the absence of zero values in the time series. In general, NCEP variables reproduce pretty well the surface observations; the GCM model tends to overestimate the precipitation, especially in November-December and spring. Finally, after verification of the GCM model, the CanESM2 data were downscaled using the SDSM for precipitation and temperature under the RCP scenarios in future period 2017-2100.
The Mann-Kendall trend test was performed for simulated and observational data to investigate the trend in climate parameters in future period compared with the base period (Table 2). Based on the results of the Mann-Kendall test, observational and modelled precipitation data under two scenarios, RCP2.6 and RCP8.5, showed no significant trend, while the data generated by the RCP4.5 scenario showed a significant decreasing trend. The results of the Mann-Kendall's temperature test showed a significant, increasing trend for temperature (observation and modelled data under all three scenarios). Also, the percentage changes in annual average modelled series under the RCP scenarios as compared with the observation values showed that the annual averaged precipitation according to the RCP2.6 and RCP 8.5 scenarios would be 62.28 and 69.00 mm, which means it would increase, respectively, by 9.18% and 17.2% compared with the observation period. In contrast, the annual average precipitation under the RCP4.5 would be 58.39 mm, which is about the same as in the observation period. The average annual temperature according to the RCP2.6, RCP4.5 and RCP8.5 scenarios would be 20.59, 21.03, and 21.22 C with percentage increases of 4.5%, 6.7% and 12.8%, respectively. Figure 5 shows the error bar of the monthly average precipitations. The average precipitation from December to March (from January to March) decreased as compared with the historical period based on the RCP4.5 and RCP2.6 (RCP8.5) scenarios, while it increased in other months. The results of simulation of precipitation based on the RCP scenarios are consistent with the results of Dastorani et al. (2011), whose study was performed for the A1 scenario. Figure 6 shows the error bar of the monthly average temperature. Based on all scenarios, the average temperature in future will increase compared with the historical period, the change being higher in the warmest months, May-August, and lower from October to December.

| Calculation of the SPI and SPEI time series in the historical and future periods
Gamma and Generalized Pareto distributions were fitted to the data for the calculation of the SPI and SPEI at six, nine and 12 months in the present and future periods, respectively. Drought characteristics, such as severity, duration and frequency, were determined according to run theory by choosing the threshold level −1 for the SPI and SPEI time series. The negative values of SPI/SPEI represent a period where the monthly precipitation is less than its long-term average, while positive values of SPI/SPEI indicate that the monthly precipitation is greater than the longterm average. According to run theory, 15-73 drought events were obtained based on the SPI and SPEI in the historical and future periods, respectively. Figure 7 shows that, with increasing SPI time scale, the frequency of the dry period decreases, which is accompanied by an increase in the severity and duration of drought. According to the SPI6, the most severe drought had a severity of −7.94 and a duration of six months, and the SPI12 had a severity of −22.08 and a duration of 11 months (Table 3). The inclusion of average temperature for calculating dry periods makes the SPEI > SPI, so the former can represent better drought characteristics (Figure 8). The most severe drought based on SPEI6 had a severity of −30.43 and a duration of 14 months. Also, the most severe droughts based on the SPEI9 and SPEI12 had a severities of −78.54 and −130.94 and durations of 44 and 70 months, respectively (Table 4). According to the SPI, the study area will experience a more severe and longer drought in future according to the RCP scenarios (Table 3). The percentage increase of the severity of drought generated by the RCP2.6, RCP4.5 and RCP8.5 scenarios according to the SPI6 was 25.8%, 42.4% and 56.4%, respectively, as compared with the historical period. The increase in the SPI9 and SPI12 was 28.6%, 28.0% and 47.1% and also 47.8%, 7.2% and 58% for the three scenarios, F I G U R E 5 Error bar (a) and shaded error (b) of monthly average precipitation in the historical and future periods F I G U R E 6 Error bar (a) and shaded error (b) of monthly average temperature in the historical and future periods respectively. Among the different scenarios, the severity and duration of the drought resulting from the RCP8.5 were higher than for the RCP2.6 and RCP4.5. The results showed that the drought characteristics derived from the SPEI, under the RCP scenarios, were different from the SPI (Table 4). In SPEI6, the drought severity generated by the RCP2.6, RCP4.5 and RCP8.5 scenarios increased by 12.34%, 21.81% and 42.46%, respectively, whereas in the SPEI9 and SPEI12 under the RCP scenarios, the droughts in future were less intense than in the observation period. Similar to the SPI, the drought derived from the SPEI under the RCP8.5 scenario had a higher severity and longer duration than under the other scenarios.

| Bivariate drought analysis in the present and future periods based on the SPI and SPEI indices
The purpose of this section is the determination of the marginal distribution functions of drought severity and duration. Among the three time scales of six, nine and 12 months, the six month scale for both the SPI and SPEI  F I G U R E 8 Standardized precipitation evapotranspiration index (SPEI) time series for historical and future droughts at six, nine and 12 month time scales was selected for bivariate analysis because of higher drought frequency than the other time scales. The relationship between the characteristics of droughts showed a significant correlation between the severity and duration of the drought in both the present and future periods (τ > 0.8). The SPEI6 series in the historical period (τ = 0.91) and under the RCP8.5 (τ = 0.94) had the highest correlation and the SPI6 series under the RCP8.5 (τ = 0.082) had the lowest correlation. A total of 20 univariate distribution functions were fitted to drought severity and duration to determine the marginal distribution functions, and the parameters of these functions were estimated using the maximum likelihood method. The Johnson SB, Generalized Pareto, Log-Logistic, Gamma and Log Normal functions were selected as the best fit for drought severity, and the Normal, Log Pearson 3, Gamma, Generalized Pareto, Log-Logistic and Exponential functions as the best fit for drought duration based on statistical tests including Kolmogorov-Smirnov, Anderson Darling and χ 2 (Thilakarathne and Sridhar, 2017;Modaresi Rad et al., 2017) (Tables 5 and 6).

| Fitting the copula functions to drought severity and duration
To fit the copula functions to drought severity and duration, the parameters of the functions for both indicators (SPI6 and SPEI6) were estimated in the present and future periods using the parametric and nonparametric methods ( Table 7). The parameters of the copula functions express the dependence structure between drought severity and duration. Based on the results obtained from the selection criteria of the best copula function, three functions (the Frank, Gaussian and Gumbel copulas) were selected among the five fitted copulas to investigate the joint behaviour of the drought severity and duration in the present and future. Thus, the Frank copula for the SPI6 under the RCP2.6 and RCP4.5 scenarios and the SPEI6 under the RCP4.5 scenario, the Gaussian copula for the SPI6 and SPEI6 in the historical period and under the RCP2.6 scenario, and the Gumbel copula for the SPI6 and SPEI6 under the RCP8.5 scenario were selected because they had the lowest AIC, BIC and NRMSE criteria ( Table 7). The selection of three different copula functions showed that the joint behaviour of the drought series is different from each other. To investigate the relationship between severity and duration of the observed drought and those generated by copulas, the marginal distributions of the drought severity and duration were simulated by the Gaussian, Frank and Gumbel copula functions, as the best copulas, using the Monte Carlo algorithm (Acciolya and Chiyoshi, 2004) and compared with the marginal distributions of severity and duration of the observed drought. The results are shown in    (Figure 9c versus f); this feature follows the tail behaviour of the Gumbel copula (Trivedi and Zimmer, 2005). Figure 9 shows that the marginal function values of the severity and duration of the SPI6 under the RCP8.5 scenario have a higher correlation in the upper tail than the lower tail, which corresponds to the marginal distributions generated by the Gumbel copula for this series (Figure 9f). The behaviour of the SPI6 series in the historical period ( Figure 9b) and the SPEI6 under the RCP4.5 scenario (Figure 9a) is the same in both tails. In other words, the upper and lower tails have the same correlation, that is, this follows the behaviour of the Frank and Gaussian copulas, which have a symmetric dependence in both tails. The Gaussian copula, owing to its exponential behaviour, has the same correlation in both tails (Trivedi and Zimmer, 2005). Figure 9a,b shows that the marginal functions of severity and duration of the SPI6 in the historical period ( Figure 9a) and the SPEI6 under the RCP4.5 scenario ( Figure 9b) have a symmetric correlation in both lower and upper tails. This attitude of the SPI6 and SPEI6 series follows the behaviour of the marginal functions generated by the Gaussian (Figure 9d) and Frank copulas (Figure 9e). 3.5 | Bivariate analysis of the severity and duration of the present and future droughts The bivariate drought return period in the present and future periods was calculated using the Frank, Gaussian and Gumbel copulas, which can show the joint behaviour of the drought characteristics. The average interval time E(L) between droughts was calculated by taking into account the probability of the occurrence and nonoccurrence of droughts (Kendall and Dracup, 1992;Mathier et al., 1992): where p 1 = 1/number of drought events; and p 0 = 1/number of non-drought events. Based on E(L), if the number of droughts in the series is higher than that of wet periods, the average interval time between the drought events will be lower. For example, in the SPI6 series under the RCP4.5 and RCP8.5 scenarios, the numbers of drought events are 152 and 176, respectively.
According to E(L), the average interval time between the drought events is 8.6 and 7.6 months, respectively. Figure 10 shows the joint return period of severity and duration of drought in both T AND SD and T OR SD based on the SPI6 in the present and future. As compared with the present period, the drought return period T AND SD based on the SPI6 will increase under the RCP scenarios. As a consequence, the maximum drought return period in the present period is 50 years, while in future, under the RCP2.6, RCP4.5 and RCP8.5 scenarios, it will be 200, 80 and 600 years, respectively. On the other hand, the results showed that a drought event with a similar return period in future under climate change occurs with a higher severity and duration than in the historical period. For example, a drought event with a return period of 10 years corresponds in the historical period to severity of −6.11 and 5 months' duration, while in future under the RCP2.6 scenario to a severity of −9 and a duration of 5 months; under the RCP4.5 scenario to a severity of −6.9 and a period of 5 months; and under the RCP8.5 scenario to a severity of −5.8 and a period of 5 months. The drought return period T OR SD < T AND SD (Durante and Salvadori, 2010;Salvadori et al., 2011;Madadgar and Moradkhani, 2013). In other words, the joint return period of drought when both the variables, severity and duration, exceeded a certain value (T(S ≥ s and D ≥ d)) was longer than either variables, severity or duration, exceeded a certain value (T(S ≥ s or D ≥ d)). In contrast to the return period T AND SD , the drought return period T OR SD will be longer than the historical period only under the RCP8.5 scenario and will be reduced under the RCP2.6 and RCP8.5 scenarios. Therefore, the maximum return period of T OR SD in the historical period is 25 years, and in future under the RCP2.6, RCP4.5 and RCP8.5 scenarios, it will be 10, nine and 250 years, respectively. Figure 11 shows the joint drought return period T OR SD and T AND SD based on the SPEI6. According to the SPEI6, the joint return period T AND SD under the RCP scenario will be longer as compared with the historical return period. In particular, the maximum joint return is 13 years in the historical period, while in future under the RCP2, RCP4.5 and RCP8.5 scenarios it will be 200, 106 and 36 years, respectively. According to the SPEI6, a drought event with a return period of approximately 10 years in a historical period occurred with a severity of −24 and duration of 12 months, while the drought with the same return period in future under the RCP2.6 scenario has a severity of −19 and duration of 6 months, under the RCP4.5 scenario with a severity of −22 and duration of 13 months; and under the RCP8.5 scenario with a severity of −30 and duration of 17 months. Figure 11 shows that, similar to the SPI, the return period T OR SD < T AND SD . However, according to the SPEI, in contrast to the SPI, the return period of the future drought in both cases of the return period (T AND SD and T OR SD ) increases relative to the historical period. Therefore, the maximum drought return period T OR SD based on the SPEI6 in the historical period is nine years and in future under the RCP2, RCP45 and RCP8.5 scenarios will be 98, 20 and 32 years, respectively.

| Conditional return period
The drought return period based on the two characteristics, severity and duration, in addition to the two cases T AND SD and T OR SD can also be expressed as conditional (Shiau, 2006). In this case, the conditional return period (T(D | S ≥ s)) is determined for drought duration given a drought severity that exceeds a certain threshold (Equation 15): where T DjS≥s 0 is the conditional return period for D given S ≥ s. Figure 12 shows the conditional return period of drought duration based on severity larger than the certain threshold in the historical and future periods under the RCP scenarios for the SPI6 and SPEI6 series. In the present study, drought severity threshold(s) was determined based on the maximum and minimum severities. For example, the return period of a drought event according to the SPI6 in a historical period with a duration of four months and a severity > 6 (< −6) is 15 years. At this level (S ≥ 6), the return period for drought duration exceeding four months, given a drought severity > 6 under the RCP2.6 scenario, is 13 years, under the RCP4.5 scenario is 16 years, and under the RCP8.5 scenario is 23 years. Also, according to the SPEI, the historical return period for drought duration exceeding four months, given a drought severity > 6, is three years, F I G U R E 1 0 Joint return period of drought T AND SD (left) and T OR SD (right) in the historical period (a) and under representative concentration pathway scenarios (RCP2.6 (b), RCP4.5 (c), RCP8.5 (d)) based on the standardized precipitation index SPI6 under the RCP2.6 scenario is six years, under the RCP4.5 scenario is six years, and under the RCP8.5 scenario is three years. The results showed that the drought return period in the joint and conditional cases is longer according to the SPI as compared with the SPEI. In other words, according to the SPEI, the study area will experience more severe droughts and a shorter return period in future.
F I G U R E 1 1 Joint return period of drought T AND SD (left) and T OR SD (right) in the historical period (a) and under representative concentration pathway (RCP) scenarios (RCP2.6 (b), RCP4.5 (c), RCP8.5 (d)) based on standardized precipitation evapotranspiration index SPEI6

| CONCLUSIONS
Most of Iran has a dry and semi-arid climate, and it is very much affected by climate change. A decrease or increase in climate parameters, such as precipitation and temperature, is affecting extreme events such as droughts. In the present study, the meteorological drought in Yazd province, located in the centre of Iran, as an example of a dry area under water stress and long droughts (Dastorani et al., 2011) was studied in the present and future periods under climate change conditions. The standardized precipitation index (SPI) and standardized precipitation evapotranspiration index (SPEI), considered as the most used meteorological drought indices, were employed to study the drought condition in future under climate change. Using copula theory, the joint behaviour of drought duration and severity was analysed in terms of the SPI and SPEI. Based on the selection criteria of the best copula, three copula functions were selected, that is, Frank, Gaussian and Gumbel, for the analysis of the joint behaviour of drought characteristics. The six month time scale was selected for bivariate analysis because of the higher frequency of droughts than other scales in both indices, the SPI and SPEI. The SPEI, which considers the effect of evaporation on drought, because of including the temperature parameter in addition to precipitation for calculating drought periods, showed more severe droughts with higher duration in comparison with the SPI (Homdee et al., 2016). In fact, the use of the temperature parameter allows the consideration of the impact of climate change on the drought situation in the area, and the drought is determined better than the SPI in terms of global warming. According to the results of both indices (SPI6 and SEPI6), the study area will experience a more severe and longer drought in future; in particular, based on the RCP8.5 scenario, more severe droughts will occur in future. Also, according to the SPEI6, the severity of droughts will increase as compared with the SPI6 (Homdee et al., 2016). The joint and conditional return periods of the drought characteristics were calculated using the Gaussian, Frank and Gumbel copula functions for the SPI6 and SPEI6 series in the historical and future periods. The results showed that the joint return period of drought when the severity and duration exceeded a certain value (T AND SD ) was longer than that either variables, severity and duration, and exceeded separately a certain value (T OR SD ). According to the SPI6 and SPEI6, the study area will experience droughts with a longer return period in future as compared with the historical period. The joint and conditional return periods are a useful tool for managers and water resource planners to provide useful information for risk assessment (Shiau, 2006;Song and Singh, 2010). By knowing the joint return period of drought in future, one can improve the policies of water resources management since current policies may not be able to reduce the negative effects of droughts. In general, the results have important implications for the management of environment risk in a long period.