The Impact of Climate Change on Hydro-Meteorological Droughts Using Copula Functions

Climate change has made many alterations to the climate of earth, including hydro-climatic extreme events. To investigate the impact of climate change on hydro-meteorological droughts in the Kamal-Saleh dam basin in Markazi province, Iran, proportional to future climate conditions, a new and comprehensive index was developed with the aim of accurately estimating drought in a more realistic condition. This aggregate drought index (ADI) represented the main meteorological and hydrological characteristics of drought. Temperature and precipitation projections for future climates were simulated by five CMIP5 models and downscaled over the study area during 2050s (2040–2069) and 2080s (2070–2099) relative to the baseline period (1976–2005). By fitting five univariate distribution functions on drought severity and duration, proper marginal distributions were selected. The joint distribution of drought severity and duration was chosen from five types of copula functions. The results revealed that in future, severe droughts are expected to frequently occur in a shorter period.


Introduction
Today, the unprecedented environmental and socio-economic impacts of climate change have become undeniable. Recent studies have suggested that changes in the frequency of extreme events may be accompanied by climate change. In general, an increase in droughts, especially in hot seasons at lower and middle altitudes, where there is the highest water demand, is expected in the near future. To monitor droughts, many studies have defined them as "lack of rainfall" based on precipitation data (e.g., Jain et al. 2015;Kundu et al. 2020).
It is now evident to experts that a univariate method is not deemed to be the best method for describing characteristics of draughts (Hao and Singh 2015;Kwon et al. 2019). Droughts constitute multivariate events because any fluctuations in climate variables can increase or decrease their severity (Kundu et al. 2020). Moreover, over the years, efforts have been made to develop multivariate indices based on a combination of different drought indices to provide insight into the conditions and characteristics of this climatic event (Beersma and Buishand 2005;Li et al. 2015). To quantify multivariate indices, various methods and theories, such as establishing combinations of precipitation and evapotranspiration variables in different statistical models (Xia et al. 2014), copula functions (Nelsen 2006;Salvadori and Michele 2007;Vergni et al. 2015;Mortuza et al. 2019;Tahroudi et al. 2020), scalogram model (Sough and Mosaedi 2013;Sough et al. 2018), Principal Component Analysis (PCA) (Keyantash and Dracup 2004;Rajsekhar et al. 2014), entropy theory (Hao and Singh 2015), and Markov chain method (Rezaeianzadeh et al. 2016) have been developed. The PCA is a feature extraction and information compression method with a long history in linear algebra and vector analysis (Sough et al. 2018). Through mathematically formulating the approach, this method is completely pliable to develop an aggregate drought index (ADI) from a combination of effective parameters on meteorological, agricultural, and hydrological indices to comprehensively monitor droughts. Studies have also reported that the ADI structure provides a clear and purposeful way to describe the severity of droughts and can be easily used to express their characteristics in the actual projects (Bazrafshan et al. 2014;Li et al. 2015).
Given that droughts are random and multivariate events, their characteristics (severity and duration) are interdependent, meaning that these characteristics do not change independently, and any changes in one variable influence the other one. Furthermore, the bivariate copula function between these two main characteristics reveal a significant correlation between them so that by employing this function, the return period of drought events are determined (Yusof et al. 2013). Kwak et al. (2015) tried to study the impact of climate change on future hydrological droughts. They applied the copula method to analyze the joint probability distribution of droughts. Azam et al. (2018) also estimated the values of four identified homogeneous regions by combining important drought variables (duration and severity) based on the standard precipitation index (SPI). The regional frequency analysis of droughts was conducted by evaluating the probability distributions and copulas. The Pearson type III and Kappa marginal distributions simulate the regional drought variables better. Bivariate stochastic simulation of the selected copulas further illustrated that the behavior of simulated data may change when the correlation between drought characteristics is considered so the joint distributions were used to estimate conditional probabilities (Mirabbasi et al. 2012).
In the present study, drought conditions due to climate change in the Kamal-Saleh dam basin in Markazi province, Iran, was projected by adopting a new aggregate drought index (ADI) calculated by PCA method from two different types of meteorological and hydrological drought index simultaneously for the present and future periods. Previous research has mostly either focused on the drought index separately or utilized different multi-index without using copula function for drought analysis. Therefore, what has differentiated this study from prior related ones are developing a hydro-meteorological index with minimum climatic parameters and estimating the joint probability distribution as well as return periods of drought severity and duration by applying the bivariate copula method to successfully monitor future droughts in the region under consideration.

3 2 Study Area and Data
This research was conducted on Kamal-Saleh dam basin, one of the most important water supplies in Markazi province, Iran (Fig. 1). It is located in the southwest of Markazi province and northeastern of Lorestan province, between 49º 04 ̍ 02 ̍ ̍ to 49º 27 ̍ 11 ̍ ̍ E and 33º 33 ̍ 13 ̍ ̍ and 33º 33 ̍ 55 ̍ N, with a catchment area of 699 km 2 . The maximum, minimum, and average height of the basin is 2960 m, 1840 m, and 9490 m, respectively. Since Markazi province is located in the arid and semi-arid regions of Iran, Fig. 1 The study basin and stations The Impact of Climate Change on Hydro-Meteorological Droughts… 3971 therefore it has always been at the risk of drought and prone to the perceptible impact of climate change.
In the present research, the required data, including daily temperature, precipitation, and monthly streamflow from 1976 to 2005 were retrieved from the Kamal-Saleh basin stations, shown in Fig. 1. Using GCMs under RCP scenarios, climate variables, such as temperature and precipitation were simulated in future periods of 2050s (2040-2069) and 2080s (2070-2099). The RCP scenarios are consistent with a wide range of possible changes in future anthropogenic (i.e., human) greenhouse gas (GHG) emissions, considering the atmospheric concentrations. On this account, two types of RCP scenarios were adopted in this study: • RCP4.5: stabilization without overshoot pathway to 4.5 W/m 2 at stabilization after 2100. • RCP8.5: rising radiative forcing pathway leading to 8.5 W/m 2 in 2100. (IPCC 2014) The obtained results were downscaled by SDSM for the study basin. Moreover, to investigate the impact of climate change on the inflow to the Kamal-Saleh dam reservoir, the IHACRES rainfall-runoff model (a monthly conceptual model) was applied to provide a framework with minimum required parameters, which reduces the uncertainties inherent in parameter estimation.

Aggregate Drought Index (ADI), Standardized Precipitation Index (SPI), and Stream Flow Drought Index (SDI)
IDA's approach provides a clear and objective approach to describing drought intensity. This index, which can adequately reflect the behavior of hydro-meteorological drought, is recommended as an integrated index for assessing and tracking droughts. Keyantash and Dracup (2004) proposed the flowchart for processing the index on a monthly scale. For this purpose, the monthly values of SPI and SDI are first calculated by Eqs. (1) and (2) for each time series. Then, the values of the indices are adjusted for each month in different years in the form of a matrix with specified dimensions (number of months in the two indices used). On the other hand, the SPI was developed by McKee et al. (1993) and many drought researchers have acknowledged the flexibility of this index as well as its ability to monitor different aspects of drought. According to this index, the drought occurs when the SPI takes negative values and ends when the SPI becomes positive. The SPI is calculated as follows: where n is the number of months for which cumulative rainfall is calculated; P 0 is normalized precipitation for the current month; P −i is normalized precipitation for the last month; μ n is average cumulative precipitation for n months, and δ n is standard deviation of precipitation for n months. (1) The SDI calculated from river discharge data can provide a good estimate of hydrological droughts. This index can be expressed in the following manner (Nalbantis and Tsakiris 2009): where V i,k is runoff value for i th month and k th scale; V k is average runoff for k th scale; S k is deviation for k th scale.
One of the most common tools to reduce a data set involving a large number of variables to a new smaller data set is PCA (Wilks 2011). According to this multivariate statistical technique, the eigenvectors were first computed for each set of indices arranged in a matrix. The principal component values were then estimated from a linear combination of the input indices for each month and the values of the ADI were calculated. The PCA can be also employed to distill the critical information from different kinds of available drought-related variables or indices. According to this concept, the ADI as a multivariate drought index is formed to include all physical forms of drought, such as hydrological and meteorological (Keyantash and Dracup 2004;Hao and Singh 2015). By the PCA method, the ADI is the first principal component that is normalized via monthly time series of standard deviation. Since the ADI is a standardized index, it follows the same drought classification criterion as SPI. This way, a drought event continues as long as the ADI is continuously negative and becomes -1 or less and ends when ADI is positive. It should be noted that for some drought events, although the SPI value has never been lower than -1, due to the long duration of these events, their cumulative amount and effect is greater than other droughts whose SPI is less than -1.
Generally, acute problems of water resources during the drought period are related to such long-term events. In this study, according to the recommendation of Loukas and Vasiliades (2004), the drought event was defined as a period that the SPI values, and consequently the ADI values are less than zero. Therefore, positive values of ADI indicate overall wet conditions, while negative and zero values suggest overall dry and normal conditions, respectively.
In the present study, furthermore, Run theory was used to define drought characteristics. Based on this theory illustrated in Fig (2) The Impact of Climate Change on Hydro-Meteorological Droughts… 3973 In addition, as shown in Fig. 2, the drought duration (D) is defined as the duration of continuous drought.

Selection of Appropriate Marginal Distribution
The classical study by Sklar (1959) had the imperative results pertaining to this topic by introducing the notion and the name of a copula, and proving the theorem bearing the same name now. From another perspective, copulas are functions that combine or "couple" multivariate distribution functions with their one-dimensional marginal distribution functions. On the other hand, copulas are multivariate distribution functions with one-dimensional margins consistent over the range (0, 1) (Nelsen 2006).
The present study used copula functions to create a severity-duration bivariate distribution. So, technically, the copula approach divides the problem into two mutually independent steps concerning fitting the marginal distributions and finding the dependence structure among the margins (Vo et al. 2020). In this study, five univariate distribution functions were fitted to the severity and drought observation data using the Kolmogorov-Smirnov test. The distribution functions adopted in this research are listed in Table 1. The maximum likelihood estimation (MLE) was also used to estimate the parameters of f (x) (parameters of marginal distributions), and the Kolmogorov-Smirnov (K-S) test method was employed to evaluate the fitting degree. Proper copula functions were also used to integrate both D and S.

Copula Functions
In this study, five bivariate copula functions, including Clayton, Frank, Gumbel, Joe, and Guassian were considered to fit the marginal distributions and create the joint probability of drought severity and duration. Moreover, estimating the parameter of the copula functions was needed to fit these functions to drought duration and severity. There are several methods to estimate the joint dependence parameter in the data. In this study, the parameter was estimated by the log-likelihood function. Further, the best-fitted Copula was selected based on maximizing the log-likelihood function of the copula model (Nelsen 2006). Copula families adopted in this research are listed in Table 2.

Joint Probabilities of Droughts
The probability that both severity and drought characteristics exceed a certain threshold provides useful information to improve the water resources management under drought conditions. This probability cannot be calculated through a univariate analysis of severity and duration of drought, whereas it can be easily obtained by the copula functions (Shiau 2006). The joint distribution function of drought duration and severity is defined by Copula function as follows: where P (D ≥ d, S ≥ s) is the probability that both severity (S) and drought duration (D) exceed a certain threshold;F d (D) and F s (s) are the cumulative drought duration and severity distribution functions, respectively, and C F d (D), F s (s) introduces the joint probability calculated by copula function. Duration and severity are two characteristics that describe droughts as bivariate events.

Correlation Analysis of Variables
There are several methods to investigate the correlation between drought characteristics. The dependence measures can be assessed using Pearson's linear correlation r, Spearman's rank correlation ρ, and Kendall's τ (for more details, refer to Chen and Guo 2019).

Bivariate Drought Return Period
Planning and management of water resource systems under drought conditions requires estimates based on the return period of drought. Shiau and Shen (2001) theoretically derived the return period for droughts with severity greater than or equal to a given value, sʹ, according to the predicted drought inter-arrival time and the distribution function of the cumulative severity of drought, expressed as follows: The Impact of Climate Change on Hydro-Meteorological Droughts… 3975 where L is the drought inter-arrival time; E(L) is the expected drought inter-arrival time, which can be estimated from the observed droughts, and T s is the return period defined solely by the drought severity. Similarly, the return period for drought duration greater than or equal to a certain value, dʹ, can be written as follows: where T D is the return period defined only by the drought duration. Droughts are regarded as bivariate events characterized by their duration and severity. Shiau (2006) suggested a methodology that classifies the return periods of bivariate distributed hydrologic events as joint and conditional return periods. According to the joint probability distribution of droughts, return periods for D ≥ d and S ≥ s as T DS and return period for D ≥ d or S ≥ s as T ′ DS are defined for base and future periods where F D (d) and F S (s) are the cumulative distribution functions of drought duration and severity, C is the copula function and E(L) is the average inter-arrival time.

Conditional Return Period
Similarly, the conditional return period can be prepared for two cases: (a) T D|S≥s ′ (b) T S|D≥d ′ These conditional return periods are obtained as follows (Yusof et al. 2013): where T D|S≥s ′ is conditional return period for drought duration, D, while S ≥ s ′ and T S|D≥d ′ is conditional return period for drought severity, S, while D ≥ d ′ , parameters d ′ and s ′ are specific values of duration and severity, respectively.

Selecting the Best GCM for Kamal-Saleh Dam Basin
In this study, a suitable GCM for the studied region was selected by hindcast approach whereby the performance of GCMs were evaluated through comparison of historical model simulations with observations (Farjad et al. 2019). The five most widely used GCMs were initially selected and evaluated by the hindcast approach (Table 3 and . 3). Also, to assess the performance of GCMs, the coefficient of correlation (R), root mean square error coefficient (RMSE), and mean absolute error (MAE) were applied. Based on the model comparison approach, CanESM2 was determined as one of the best performance models compared to other GCMs in the studied region.

Calculation of Drought Indices
The hydrological index (SDI) and meteorological index (SPI) are calculated on the monthly scale. The values of the indices are shown in Fig. 4 for the base period. As this figure illustrates, although these two indices are generally consistent, there are differences between them over several time intervals. One of the reasons for this difference could be the fundamental role of temperature on the rate of evaporation and transpiration in this semi-arid region. This is probably because a sizable portion of rainfall is lost as water vapor to the atmosphere by evaporation and transpiration, which greatly impact the rate of runoff and subsequently lead to hydrological drought. Another reason is likely due to the abnormally heavy rainfall in a very short time, while the rest of the month is dry (SPI > 0 and SDI < 0). ADI was calculated by the PCA method described in Sect. 3.1. The following steps summarize the calculation of this index for the base period.
Initially, the matrix X was created. In this matrix, the rows indicate the months in the base period and the columns show associated SPI and SDI values: To compute the main component of the PCA requires the production of a symmetric square matrix (R) to describe the relationship between the main data. R is the correlation 1 3 matrix between SPI and SDI Eq. (11). Here, N is the number of months in the 30-year period (360 months) and X T is the transpose of the matrix.
eigenvector array for this time series is determined as: Subsequently, ADI is defined as Eq. (13) where σ is the principal component's standard deviation. Figure 5 shows the ADI values in the base period in the studied region. According to this figure, the historical values of the ADI vary from -3.5 to 3.5. By Fig. 5, it can also be understood that the most severe drought had happened in January 1995 and October 2004, the longest duration were 11, 10 and 7 months, which had occurred from May to March 1976, January to October 1980, April to October 1992, and March to September 1993, respectively. Figures 6 and 7 show ADI values for the two future periods (2050s and 2080s) under two RCP scenarios. According to Fig. 6, drought index values reveal more fluctuation in the first 30-year period, and this period would have droughts with longer durations relative to another period. Although in the second 30-year period, the ADI values are in a larger interval, the peak points of drought are greater. Therefore, more short-term severe droughts will be seen in the far future compare to the near future. Despite the result of Figs. 6 and 7 shows that the first period would deal with more long-term severe and longer drought in comparison with the second period in the future under RCP8.5, but like RCP4.5, peak points in the second period show higher values.

Investigation of the Correlation between Drought Duration and Severity
The hydro-meteorological drought severity and duration data were extracted from Figs. 6 and 7 based on the runs theory as shown in Fig. 2. The summary description of drought events and characteristics of each scenario and time series are prepared in Table 4.
In Table 5, computed correlation metrics such as, Pearson, Kendall's tau, and Spearman's rho are prepared. Pearson's correlation coefficient only considers the linear correlation The monthly ADI (the base period) The Impact of Climate Change on  between the data and in the presence of severe dependencies, this coefficient may not be computable. Non-parametric correlation coefficients such as Spearman and Kendall Tau coefficients may overcome the defect of the Pearson coefficient (Salvadori and Michele 2007). Before establishing the joint distribution function, these three rank correlation coefficients were adopted to test the correlation between drought duration and severity. The results showed a positive correlation between the drought severity and duration variables based on three methods that are used. 1 3

Estimation of Marginal Distributions for Drought Duration and Severity
Following the recommendations of Ayantobo et al. (2017), in the selection of marginal distributions, the drought duration (D) and severity (S) were fitted with their corresponding marginal distributions. The results are illustrated in Figs. 8 and 9. In these two figures, it can be seen that the Lognormal distribution fits the drought duration and severity for all of the time series and scenarios. The parameters of these distributions were estimated by the maximum likelihood method for the base and the two future periods ( Table 6). The results showed that the Lognormal well-fitted the drought duration and severity for all-time series and scenarios. Therefore, Lognormal function, whose parameters are shown in Table 7, was selected as the optimal distribution function.

Analysis of Joint Distribution Functions of Drought Duration and Severity
After fitting the optimal marginal distribution function to the drought characteristics, five common copula functions (Table 2) were considered to analyze joint distribution functions of drought duration and severity. Moreover, to establish an optimal copula function, copula dependence parameters needed to be calculated. To this end, the maximum likelihood method was used. The AIC and BIC values are also performance criteria whose lower values imply better fitting. The most optimal copula function in each period and scenario and their parameters are presented in Table 8.

3
The selected optimal marginal distributions and optimal copula produced the copulabased joint distribution of drought duration and severity. The joint probabilities of drought severity and duration can be seen in Fig. 10. Expectedly, the contours of 0.1 to 0.9 in all cases are different. Figure 10 enable understanding of probabilities of different cases of drought severity and duration. For example, in Fig. 10b (2050s under RCP4.5), the plotted contours denotes that the probability of severe drought is less than others. Contrarily, in Fig. 10c (2050s under RCP4.5), the probability of occurrence severe and long drought is considerably more than other cases in this study. Also, for both future periods, the probability of severe drought is more than the base period under RCP8.5.

Bivariate Return Period of Droughts
The results of return periods estimated from Eqs. (7) and (8) are shown in Figs. 11 and 12. When a drought occurs with a certain return period, shorter droughts whose severity is greater than the base period are likely to occur in all cases when either duration or severity exceed certain levels. This indicates that droughts will occur with a shorter duration and greater severity in the future, especially this increase is more considerable in the second period and the RCP8.5. As a result, the maximum drought return period for D ≥ d or S ≥ s (T � DS ) during the base period is 7 years, while in future in the years 2050s under RCP4.5 and RCP8.5 scenarios, will be 22 and 17 years, respectively, and during 2080s, under both scenarios will be 3 years (Fig. 11).
Also, according to Fig. 12, maximum return period for D ≥ d and S ≥ s (T DS ) during base period is 62, in the years 2050s under RCP4.5 and RCP8.5 scenarios, will be 4000 and 540 years, respectively, and in the years 2080s, under RCP4.5, it will be 800 and under RCP8.5 will be 57 years (Fig. 12).
On the other hand, Figs. 11 and 12 illustrate the estimated return period for any given duration and/or severity exceeding a certain level. Under RCP4.5, compared with the base period, the return period will be greater, while for RCP8.5, with a certain return period, droughts will occur with a shorter duration and greater severity in the future. Furthermore, When both variables (severity and duration) exceeded a certain value, the return period was more extended; for example, in the historical period, a drought event with a severity of 5 and a duration of 7 months, will occur with a return period of approximately 14 years, while the drought with the same severity and duration in the period of 2050s will have a return period of 172 and 11 years under RCP4.5 and RCP8.5 scenarios, respectively. The corresponding values for the period 2080s will be 82 and 13, respectively. Totally, the potential droughts in the future would become more frequent and more severe, except for

Conditional Return Period
The conditional return period of drought duration and severity were defined according to Eqs. (9) and (10) for the base period and two future periods under RCP 4.5 and RCP 8.5 scenarios and the results are demonstrated in Figs. 13 and 14. Conditional return periods can provide valuable information to water resources managers for risk assessment of water resources. It is always necessary to measure the ability of water supply system in every region to meet and manage its water needs in certain drought conditions. Figure 13 shows the conditional return period of drought duration for a given drought severity (s) where s exceeds the threshold level (sʹ), for instance sʹ = 1, 5 and 10. For example, if the severity of the drought exceeds an index of 5 over a drought period of more than 3 months, the return period is approximately 65 years in the historical period, whereas over the 2050 period, the return period would be about 120 and 40 years under scenarios RCP4.5 and RCP8.5 respectively. These corresponding values are estimated to be about 90 and 30 years, respectively, during the 2080s.
In another case, the study of the conditional return period of drought severity illustrated in Fig. 14 found that during the historical period, when the drought duration exceeds 5 months and its severity is more than 3, the return period is approximately 30 years while the corresponding values would be about 40 and 37 years in the 2050s and about 22 and 20 years in the 2080s under scenarios RCP4.5 and RCP8.5, respectively. Although the plots with various sʹ values had similar patterns, but, the results indicated that for a given drought duration, the conditional return period of drought occurrences increases with increasing sʹvalues (Fig. 13). Figure 14 shows similar results. The amount of these increases will be more significant in the future periods. Also the results indicate that in future periods, the conditional return periods of a particular duration with a given severity will decrease more in 2080s compare to 2050s; in addition, those decreases in RCP8.5 compare to RCP4.5. The Impact of Climate Change on Hydro-Meteorological Droughts… 3989

Conclusion
Drought analysis was conducted for the Kamal-Saleh basin as a semi-arid area threatened by droughts. The study investigated the impact of climate change on future hydro-meteorological droughts in the study area. By using the hindcasting approach, between five mostly used GCMs, which were the most widely used in the region, CanESM2 model was chosen as the optimum model to predict the future climate change compared to other GCMs in the region by using the hindcasting approach. This model was used under RCP's scenarios (RCP4.5 and RCP8.5) to analyze drought for two future periods (2050s and 2080s). Moreover, for drought assessing and monitoring, the monthly values of SDI and SPI were calculated for the studied basin. Subsequently, the historical and future values of the ADI were estimated by the PCA method.
After examining and proving the existence of a positive correlation between the drought severity and duration, these variables were fitted with their corresponding marginal distributions and the Lognormal was selected as the best optimal distribution function.
In the next step, in order to analyze the joint probabilities and return period of drought duration and severity, among five common copula functions, the most optimal one was applied in each period and scenario to produce the joint and conditional return periods of the in the historical and future periods. In the future scenarios, the probability of occurrence more severe and longer drought was obtained significantly higher than in the baseline periods.
Additionally, according to ADI, the study area will experience more frequent drought when either duration or severity exceeding certain levels (T � DS ) for all cases, especially during 2080s under RCP8.5. But return periods estimated by duration and severity exceeding certain levels (T DS ) will significantly increase in future periods under RCP4.5, while with the same return periods, severe droughts will occur with shorter durations under RCP8.5.
The results of the conditional return period of drought characteristics also indicated that in future periods, the conditional return periods of a particular duration/severity with a given severity/duration will decrease more in 2080s compare to 2050s also under RCP8.5 compare to RCP4.5. The conditional return periods of drought severity and duration provide useful information to predict and prevent the destructive effects of drought. Each return period describes different situations; hence, depending on the type of drought danger, the type of return period should be selected.
The findings indicate that drought is a serious problem in the study area. The analysis of bivariate joint and conditional return periods of drought offers invaluable information to improve the water resources management in the region. This research also showed that in the future the studied basin will be exposed to considerable droughts as a negative impact of climate change. Therefore, it is essential to make a management plan to reduce its probable negative effects.
Author Contribution This research was conducted as a part of a master's degree thesis work of Zahra Fahimirad under the supervision of Nazanin Shahkarami. Both authors have equally contributed to the research.

Funding
The authors did not receive support from any organization for the submitted work.

Declarations
Ethics Approval Compliance with Ethical Standards.

Consent to Participate Not applicable.
Consent for Publication All Authors declare that they agreed with the content and give explicit consent for publication.

Conflicts of Interest
All Authors declare that they have no conflict of interest.