Copula-based bivariate drought severity and duration frequency analysis considering spatial–temporal variability in the Ceyhan Basin, Turkey

In recent times, copula families have been mostly employed in bivariate drought duration and severity. Archimedean copulas are mainly used conveniently among various potential families for hydrologic design and water resource management. In this study, drought events are defined by standard precipitation index (SPI3 time series) series for 24-gauge stations in Ceyhan Basin Turkey considering different elevation levels (high, average, and low). Both duration and severity are tested for their dependence with Mann–Kendall and Spearman Rho trend test at a 0.05 significant level. While Gamma and Weibull perform better than the other distributions at average elevation, Log-normal and Weibull distributions are found to fit the DD and DS series better for most stations located in low and high elevations, respectively. The best fit copula is found to construct a joint distribution function. The influence of upper tail dependencies is tested to select the most appropriate copula function. Gumbel and BB1 copulas are often regarded as the best copulas for the DD and DS series, respectively. After acquiring the best fit copula, the joint return periods are modeled for each station. Finally, considering drought risk categories (mild, moderate, severe, and extreme drought), spatial distributions of drought risk-return period are constructed. High risks are observed in the middle part of the basin based on TDS and T′DS under mild drought categories. Apart from using univariate frequency analysis, bivariate frequency analysis is performed to model the joint return period based on copula theory which is great of importance for hydrology design and water resource management.


Introduction
Drought is regarded as a complex natural reoccurring phenomenon, which may affect many areas of the world, resulting in a precipitation deficit compared with the expected or normal amount over some period (Mckee et al. 1993;Chen and Sun 2015;Esit et al. 2021). Unlike many natural disasters that occur suddenly, droughts are insidious and can produce tremendous negative impacts environmental (Mirabbasi et al. 2012;Sheffield et al. 2012), agricultural (Wang et al. 2019), and socio-economic damage (Chen et al. 2013;Tosunoglu and Can 2016). For instance, in global damages, droughts may annually cause as economic losses of about $6-$8 billion (Wilhite 2000) and collectively affected a large population (including 35% of those affected by natural disasters) than other climate-related natural disasters (She and Xia 2018).
Comprehensive drought definition has been provided by researchers (Linsley et al. 1949;Gumbel 1963;Palmer 1965;Yevjevich 1966;Wilhite and Glantz 1985;Organization (WMO) and World Meteorological Organization (WMO) 1986; Mishra and Singh 2009;Schneider 2011) due to complexity of causes and influencing factors of drought. Drought may generally be classified into four categories, i.e., meteorological, hydrological, agricultural, and socioeconomic drought (Mishra and Singh 2010). In drought analysis, a variety of indices have been developed to assess and monitor droughts. The Standardized Precipitation Index (SPI) (Mckee et al. 1993) and Palmer Drought Severity Index (PDSI) are widely used methods for drought characterization. PDSI, based on water balance accounting using precipitation and temperature parameters, is generally used within the USA. PDSI has little acceptance by many researchers due to several limitations (suspicion of water balance model, temporal scale, PDSI value) (Kao and Govindaraju 2010). Unlike PDSI, SPI, which considers the deficiency of precipitation, is based on the probability concept to classify the drier and wetter climates (Cancelliere et al. 2007). Owing to its simplicity for calculating multiple time scales, SPI is widely accepted all over the world using precipitation parameters (Guttman 1998) solely.
Drought characteristics are usually analyzed separately by univariate frequency analysis (Cancelliere and Salas 2004;Serinaldi et al. 2009). Because of complex phenomena, using one parameter cannot depict the extensive impact of drought events. Two significant randomly correlated drought characteristics are generally used in literature as severity and duration (Shiau and Modarres 2009;Shiau et al. 2012). For investigating the impact of drought, a convenient method is to employ the stochastic process and the probability theory method (Shiau 2006;Tosunoglu and Kisi 2016). However, since drought variables (e.g., severity and duration) are randomly correlated to each other, univariate analysis of drought characteristics cannot provide a significant correlation between these variables. Multivariate analysis characteristic is a better approach for assessing drought characteristics. But, multivariate models are complicated to derive joint distribution functions of drought variables (Mishra and Singh 2010). For instance, Bonaccorso et al., Kim et al., Salas et al., and Shiau and Shen (Shiau and Shen 2001;Bonaccorso et al. 2003;Kim et al. 2003;Salas et al. 2005) have suggested different theories for identifying the joint distribution of drought variables. However, due to some difficulty in mathematical derivations and obtaining their parameters by fitting the observed data, these models are not applicable (Chen et al. 2013).
Bivariate and multivariate distributions using copula functions may overcome the above-stated problems. The copula is an effective method to assess the dependency and correlation among multiple variables by linking multivariate distributions depending on their univariate marginal distribution (Sklar 1959). In recent years, copulas are widely used as a method for addressing multivariate hydrological analysis problems. For example, copulas have been employed for flood frequency analysis (Zhang and Singh 2006;Shiau 2006;Wang 2016;Alamgir et al. 2018;Durocher et al. 2018;Dong et al. 2019), rainfall frequency analysis (Zhang and Singh 2006;Kuhn et al. 2007;Kao and Govindaraju 2010), drought frequency analysis (Masud et al. 2015;Kwon and Lall 2016;Vazifehkhah et al. 2019;Wang et al. 2019;Hui-Mean et al. 2019), for predicting groundwater parameters (Bárdossy 2006), and for estimating streamflow discharge using remote sensing (Chacon-Hurtado et al. 2017). A comprehensive description of theoretical backgrounds for using copulas can be found in Salvadori et al. (2007) and Nelsen (2007).
For bivariate drought frequency analysis, Shiau (2006) and Shiau et al. (2007) constructed joint probabilities and return periods for drought severity and duration using two-dimensional copulas. Yoo et al. (2016) modeled confidence intervals of the bivariate drought frequency curve after analyzing the model performance of several copula functions. Applying the Copula-GARCH rainfall generation model, 100 realizations of 100-year-long monthly precipitation were produced. Ahmadi et al. (2018) applied a copula for estimating bivariate frequency analysis of low flow in the Dez River Basin, Iran. Song and Singh (2010a) constructed the joint probability distribution of more than two variables including drought duration, severity, and inter-arrival time, employing a trivariate Placket copula. In the other their study, they applied various metaelliptical copulas (Gumbel-Hougaard, Clayton, Frank, and Ali Mikhail-Haq) to obtain trivariate joint distributions of drought variables (Song and Singh 2010b).
A few studies have been performed for bivariate frequency analysis of any specific region in Turkey using copula functions (Vazifehkhah et al. 2019;Afshar et al. 2020;Tosunoğlu et al. 2020;Avsaroglu and Gumus 2022). There is no research applied to the Ceyhan Basin, Turkey. The novelty of the study is to determine the most appropriate copula function to obtain a joint return period of drought considering different elevation levels (high, average, and low). Evaluating copula functions according to the elevation levels of the basin allows a better interpretation of the study results. In addition, based on its agricultural, irrigation, and hydropower facilities, the Ceyhan Basin is one of the most important crucial basins for Turkey. For this reason, this paper's goal is to investigate the univariate and bivariate distributions and derive the joint probability distribution of meteorological drought in the Ceyhan Basin, Turkey. To achieve this purpose, (1) monthly precipitation data from 24 meteorological stations located in the Ceyhan Basin are gathered and the Mann-Kendall trend test is employed to determine the homogeneous precipitation basin; (2) the standard precipitation index (SPI) method is used to identify drought variable (severity and duration) and to test their stationary and randomness, Mann-Kendall and Spearman Rho Tests are performed; (3) the best fit marginal distributions of drought severity and duration are determined, respectively; (4) 10 types of copulas (i.e., Gaussian copula, Student t copula (t-copula), Clayton copula, Gumbel copula, Frank copula, Joe copula, BB1 copula, BB6 copula, BB7 copula, and BB8 copula) are applied to generate two-dimensional joint distributions; and (5) consideration with upper and lower tail dependence, the best fit copula for each stationary station, significant probabilistic specifications of droughts are derived. After all, the results are presented and discussed.

Homogeneity test
The Wald-Wolfowitz test, often known as the Runs test, assesses the data's randomness and examines whether one observation influences subsequent time series data that are cut from a specific level, which could be the mean, medium, or mode to evaluate if each value in the series is lower or higher than this level. The run number is the number of passes from one data to the next above or below a specific level. If the run is under or over long periods, the run number is small. These series might not be homogeneous (Wald and Wolfowitz 1942). The test's outcome is represented by z, the quantity of data N, the quantity of runs R, the number of values under the medium level Na, and the number of values over the medium level Nu: Wallis and Moore phase frequency test (Wallis and Moore 1941) is performed to detect deviations of time series due to randomness in the values' sequence. The test is based on sign differences (− or +), and the number of phases is determined without taking into account the first and last phases, also known as the sequence of signs. If n ≥ 10 and continuity correction is performed, a fairly good test can be based on the hypothesis that data is normally distributed; when n ≥ 25, the correction is not applied (Wallis and Moore 1941). Calculating the z-test statistic is as follows: where h is the number of phases, but only the first and last phases are not included. Normally distributed describes the z-statistic. For n ≤ 30, a − 0.5 continuity correction is added to the denominator. (1)

Definition of drought events
In this study, the standard precipitation index (SPI) developed by Mckee et al. (1993) was carried out for identifying drought events (Fig. 1). The main advantage of using SPI is statistical consistency. It also enables us to define both short-term and long-term drought effects at different time scales of precipitation anomalies. In statistics, the SPI value is equivalent to the Z-score. However, due to precipitation observations with a time scale of 12 month or less, the distribution of observation is generally considered skewed. The more appropriate distribution for precipitation data is found as gamma distribution by Thom (1958). The probability density function can be calculated for gamma distribution g (x): where x is the precipitation observation and > 0 and > 0 are the shape and scale parameters. The gamma function Γ( ) is defined as Edwards and McKee (1997) applied the maximum likelihood method for predicting and parameters: where n is the number of precipitation observations and x is the mean of x. The cumulative distribution may be defined as where q is the probability of zero and G(x) is the cumulative distribution for a desired month and time scale. Finally, the cumulative distribution H(x) is then converted to the standard normal random variable Z (mean 0 and variance 1), which represents the value of SPI (Tsakiris and Pangalou 2009). Figure 1 explains the time series of SPI and the drought process. Mckee et al. (1993) categorized droughts as a period length in which the SPI is below − 1 and SPI classifications for wet and dry events are listed in Table 1. Drought duration (D) is expressed as the number of consecutive intervals (months) where the SPI value fell below zero and drought severity (S) is defined as cumulative SPI value during drought duration, stated using the following equation as (Mckee et al. 1993) where S is the drought severity and D is the drought duration. Interval time E(L) is the time between the initiation of the drought to the beginning of the next drought (Mishra and Singh 2010).

Measure of dependence drought characteristics
In order to employ copula function, stationarity or independence of drought variables should be defined (Daneshkhah et al. 2016). In this study, Mann-Kendall and Spearman Rho, which are non-parametric trend tests that mostly used measures of dependence, were applied to calculate the serial stationarity structure of drought duration and severity. According to two test results, computed Z values were lower than the critical of ± 1.96 at a 0.05 significance level (Table 2). That is, for fitting the copula function, drought variables were displayed to be appropriate. Consequently, after proving the stationarity of random variables, the bivariate copula function can be applied. However, non-stationary distribution models would be assessed, if the data were not stationary. The Mann-Kendall's non-parametric rank correlation coefficient was also employed to drought variables. Highly correlated between variables are reasonable for applying joint distribution.

Distributions of univariate drought variables
Before constructing bivariate distribution of drought variables (duration and severity), a univariate distribution of drought characteristics must be first described. Some previous studies showed exponential for drought duration and gamma distribution for drought severity commonly used in drought analysis (Shiau 2006;Shiau and Modarres 2009;Lee et al. 2013). However, Yusof et al. (2013) have suggested that two mostly used probability distributions are not appropriate in every case well. Therefore, several distributions are used to fit the drought severity and drought duration in this study. The distributions are Log-normal, Logistic, Gamma, Exponential, Weibull, and Normal distributions. The performance of used distributions for each drought characteristic is tested by the Anderson-Darling (AD) (Stephens 1974), Kolmogorov-Smirnov (K-S) (Smirnov 1948), Cramers-von Mises (CvM), Akaike's Information Criterion (AIC) (AKAIKE 1976), Bayesian Information Criterion (BIC) (Stone 1979), and Maximum likelihood methods.

Copula functions
The copula is a powerful function that enables us to model a joint distribution function of different univariate random variables (Nelsen 2007), and it is based on the correlation between variables to link different marginal distributions. When compared to traditional multivariate distributions, copula provides great flexibility to select the univariate marginal distributions. The theoretical definition of a copula was first expressed by Sklar (1959). According to Sklar's theorem, if two random variables x and y are considered with their following the marginal distribution function F X (x) and F Y (y), respectively, then there exists a unique copula C, which links these two different marginals to model the joint distribution function, F X,Y (x, y). C is defined using the following equation as (Nelsen 2007) Equation 8 shows that a copula can define a multivariate distribution concerning a univariate distribution. Conversely, for any univariate distributions F X (x) and F Y (y) and any copula C, the function F X,Y (x, y) stated above is a twodimensional distribution function with marginal distributions F X (x) and F Y (y). Additionally, if F X (x) and F Y (y) are continuous, then C is completely unique (Shiau 2006). Then the joint probability density function defines the following: where f X (x) and f Y (y) are the density functions related to F X (x) and F Y (y), respectively, and c is the density function of C, defined as While many copula families exist including (1) Archimedean (Gumbel, Frank, Clayton, and Ali-Mikhail-Haq); (2) Elliptical (normal and t); (3) Extreme Value (Galambos, Gumbel, Husler-Reiss, etc.); and (4) other families (Farlie-Gumbel-Morgenstern, Plackett, BB1, BB6, BB7, and BB8), just Archimedean is convenient for hydrologic applications (Joe 1997;Abdous et al. 2005;Nelsen 2007). One parameter of Archimedean copula can be defined as (Nelsen 2007) where ∅ [0, 1] is the copula generating function and ∅(1) = 0 . The summary description of copula families used in the study is shown in Table 3.

Tail dependence assessments
The terms of upper and lower tail dependence were introduced by Sibuya (1960), which are related to dependence between extreme values in the upper and lower tails of the bivariate distribution. The lower positive values of the upper tail show that they are of poor relation or not strongly dependent (Reddy and Ganguli 2012). In hydrology, disregarding tail dependencies can cause high uncertainty in extreme quantile predictions, resulting in incorrect findings for hydrological design (Xu et al. 2010). On the other hand, frequency analysis of extreme hydrological events (such as droughts and floods) should be considered in terms of tail dependencies (Poulin et al. 2007). The coefficient of upper (lower) tail dependence u L is defined as where F X (x) and F Y (y) are the CDFs of random variables of drought duration and severity, respectively, and t is the Tail dependence just applies to choosing an appropriate copula family (not to the selection of marginal distribution). Previous studies showed that the best fit copula is Clayton without using tail dependence. However, with using tail dependence, while non-Gaussian copulas have a quite lower upper (lower) tail dependency, Gaussian copula does not have tail dependence (AghaKouchak et al. 2010). This paper is focused on upper tail dependence due to analyzing the occurrence of extreme events. That is, according to test results and literature, Ali-Mikhail- not take consideration (Nelsen 2007). The non-parametric upper tail dependence λ CFG u is given by the following:

Return period
A general approach to the management of the water resource system, hydrologic and hydraulic facilities, is to estimate return periods of drought characteristics (Shiau and Shen 2001). Mainly, drought return periods give crucial information under drought conditions. Univariate return of drought duration greater than or equal to a certain value can be calculated as where L is the expected drought interval time, which is mentioned in Sect. 2.1, and T D is the return period defined solely by drought duration. The return period of drought severity, greater than or equal to a certain value, can be obtained using the same formula defined as Joint drought duration and severity can be estimated in two cases: return period for D ≥ d and S ≥ s and joint return period for D ≥ d and S ≥ s which is stated by T DS and T′ DS , respectively, as follows: where F S (s) and F D (d) are the CDFs of univariate drought severity and duration and C is any type of copula, respectively.

Study area and data
The Ceyhan Basin is located between latitudes of 36°30′ and longitudes of 35°20′ in the eastern Mediterranean region of Turkey. It is bordered by the Asi in the south, the Seyhan Basin in the west and northwest, and the Euphrates in the east and northeast (Fig. 2) (Tanrıverdi et al. 2010). Ceyhan Basin, covering an area of about 20.670 km 2 , includes three central provinces, namely, Adana, Kahramanmaras, and Osmaniye. According to the rainfall regime of the Mediterranean climate type in the Ceyhan Basin, the rainiest season is winter (December, January, and February), and the least rainfall is in summer (June, July, and August). The characteristics of the Fig. 2 The location of the meteorological stations and annual depth of precipitation distribution in the Ceyhan Basin Mediterranean climate are generally seen in the basin. Ceyhan Basin is under the influence of different pressure centers. However, the influence of these pressure centers varies during the year. The climate of the Ceyhan Basin is an arid environment, where rainfalls are infrequent with short duration in summer. The annual depth of evaporation is low, varying from 0 to 325 mm. Evaporation increases significantly in summer, where significant evaporation occurred in June, July, and August. In contrast, it decreases in winter, especially in December, January, and February. The Ceyhan River's total length is 425 km, its annual discharge is 82.9 m 3 /s, and the basin yield is 10.7 L/h/km 3 . The total annual rainfall is seen in Kozan at least 842 mm, and the lowest rainfall area of the basin is in Elbistan, with 395.7 mm. The average annual temperature is at most in Kozan (19.3 °C) and at least in Göksun (8.9 °C) (Eris et al. 2019;Yuce et al. 2019;Yuce and Esit 2021).
A total of 24 gauging stations which were distributed evenly were selected in the Ceyhan Basin and its surrounding areas (Fig. 2). The monthly precipitation data were obtained from the National General Directorate of Meteorology. The distribution of annual depth of precipitation data is also shown in Fig. 2. According to the figure, while the highest distribution annual depth of precipitation is seen lower basin as 381 mm, the lowest distribution depth of precipitation occurs upper basin as 894 mm. The statistical parameters for the 24 rain gauge stations are also presented in Table 4. Gauging stations are selected at different branches throughout the Ceyhan Basin to get more accurate information about basin characteristics.

Results
Homogeneity test results considering WW and WM approaches are given in Table 5. All stations indicate homogeneity according to a 95% confidence interval except two stations (17,870 and D20M020) that is nonhomogeneous by the WW method. Stations 17,870 and D20M020 are located in high and low elevations in the basin, respectively.

Marginal distribution of drought variables
Firstly, drought characteristics are calculated from monthly precipitation data for each station using SPI demonstrated in Fig. 3. Stations 17,868, 17,255, and D20M014 as an example due to restricted paper are selected for high, average, and low elevations, respectively. SPI values below the red line indicate dry time, while values above it indicate a wet period. Before determining the best fit copula, the univariate marginal distributions should be determined for each drought variable. Tables 6 and 7 show the number of marginal distributions and evaluation of the fitted marginal distribution and GoF test results for the DD and DS series at a 3-month time scale. As stated above, in Sect. 2.2, some researchers found that the other distribution families such as Weibull and Log-normal can be more appropriate than the most used distribution functions. According to the test result performance, Weibull and Log-normal distributions for drought duration, and Weibull and Gamma distributions for drought severity, are more convenient, respectively. In a specific view, most stations located in low and high elevations are found to fit the DD and DS series better using Log-normal and Weibull distributions, respectively, while Gamma and Weibull distributions show better performance rather than the remaining distributions in average elevation. The parameters for each distribution are estimated using the performance of the test mentioned in Sect. 2.2. An example frequency analysis of selected appropriate marginal distributions, Weibull and Log-normal for duration and gamma for severity, are shown in Fig. 4. The probability density function (PDF) of Weibull is defined as where s | d is the drought severity | duration and η and are the scale and shape parameters, respectively. The probability density function of the gamma distribution is as follows: where Γ is the gamma function and and are the shape and scale parameters. The probability density function of Lognormal is expressed as

Estimation of joint distributions
After getting drought characteristics with fitted marginal distributions, 10 copula functions, including Gaussian, Student t, Clayton, Gumbel, Frank, Joe, BB1, BB6, BB7, and BB8 copulas, were tested to determine the best fit copula for modeling the dependence between drought duration and severity. Their performance for each 3-month SPI value was tested by Akaike's Information Criterion (AIC), Bayesian Information Criterion (BIC), and Maximum likelihood methods. Also, visual tests were generated and simulated 10,000 random pairs (u 1 , v 1 ) for the selected best copula (Fig. 5). The other main focus for selecting the best fit copula is to evaluate tail dependence. Descriptive information is stated in Section 2.6. The comparison of simulated copula functions (Fig. 5) versus observed drought duration and severity shows that higher dependence is not the lower tail between u 1 and v 1 . That is, higher dependencies indicate the upper tail. While test performance results demonstrate Ali-Mikhail-Haq u = 0, L = 0 , Frank u = 0, L = 0 , Galambos u ≠ 0, L = 0 , and Plackett u = 0, L = 0 functions are not convenient for the bivariate distribution function according to tail dependence and Kendall tau results. Table 8 illustrates that the Gaussian copula seems appropriate for selected  Table 5. Drought duration and severity are found well fitted as the Weibull distribution (parameters; = 1.47 and = 3.95) and Gamma distribution (parameters; = 0.82 and = 0.27), respectively. Univariate for each drought event is obtained using Eqs. 23 and 24. Interarrival time E(L) is estimated as 0.618 month considering drought time between the initiation of the drought to the beginning of the next drought. After this process, return periods of 10, 20, 50, 100, 200, and 500 years evaluated by drought duration and severity are summarized separately as shown in Table 10. The joint return period for drought events can be estimated  1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 Fig. 6 for different elevations. Comparison of simulated random data (gray data) and observed data (red data) using a fitted BB1 copula for station 17,255, b fitted Gumbel copula for station 17,868, and c fitted BB1 copula for sta-tion D20M014 (FDD and FSS are u 1 (CDF; drought duration) and v 1 (CDF; drought severity), respectively) Table 9 indicates the best fit copula family based on drought parameters for all used stations. According to Table 9, Gumbel and BB1 copulas are generally selected as more appropriate for DD and DS series, respectively. To be specific, while BB6 copulas can be seen in high and average elevations, Joe copula is detected only in average elevation over the basin. BB7 copula is observed in low and high elevations. In this paper, the joint return period  (Table 10). Different elevations indicate that drought events should be evaluated based on derived copula. Hydrologic facilities in terms of water resource system and the significant drought characteristics can be considered the possible maximum duration likely to happen over social and economic life (Mishra and Singh 2010). After evaluating the joint return period of drought duration and severity for each station, the spatial distributions of T DS and T′ DS for return periods under case 1 and case 2 are presented in Fig. 7. Drought risk categories are described as mild, moderate, severe, and extreme droughts considering a specific threshold for drought duration and severity, respectively (Chang et al. 2016). The percentile method (25%, 50%, 75%, and 95%), which is mentioned by Chang et al. (2016) in detail, was employed for the marginal distribution of drought duration and severity series considered four drought risk categories.
According to the SPI3-time scale scenario, compared to the southern regions, the northern regions, under all drought categories, especially D20M001 and D20M011, higher return period T DS and lower risk are observed. But, T′ DS indicates a lower return (high risk) period in the southern part of the region. In addition, the western part shows a tendency to decrease the risk from mild to the

Conclusions
Since drought is a complex phenomenon, understanding drought variable is great of importance for water resource management. In this study, bivariate drought analysis was applied based on copula functions. The previous studies mainly focused on how drought risk assessments are evaluated by using a single variable. However, drought should be assessed by employing multiple variables together due to affecting various parameters. Hence, bivariate copula functions are essential to construct a link between drought variables. According to the results, all stations show homogeneity based on WW and WM tests at a 95% significance level except stations 17,870 and D20M020. In addition, there is no trend detected in DD and DS series at 3-month time scale. Although many studies are using the copula theory, there are no studies based on the elevation (high, low, and average) of the basin. When the different elevations are taken into account, the most appropriate marginal distributions that fitted both the DD and DS parameters are more clearly understood. For drought severity, the Weibull and Gamma distributions are more convenient, whereas the Weibull and Log-normal distributions are better for drought duration. Log-normal and Weibull distributions are found to fit the DD and DS series better for most stations located in low and high elevations, respectively, while Gamma and Weibull distributions perform better than the other distributions at an average elevation. These results are appropriate with Chen et al. (2013) and Tosunoglu and Can (2016) papers. After acquiring drought duration and severity, they are fitted separately by the different probabilistic distribution functions. Ten copula functions, including Gaussian, Student t, Clayton, Gumbel, Frank, Joe, BB1, BB6, BB7, and BB8 copulas, are considered to determine the best fit copula for modeling the dependence between correlated drought duration and severity. Apart from employing several statistical tests (Akaike's Information Criterion (AIC), Bayesian Information Criterion (BIC), and Maximum likelihood methods), the performance of different copula classes is evaluated in terms of upper tail dependence. After performing statistical and tail dependence tests for modeling the best copula function for each station, based on the joint distribution function, the drought severity-duration-frequency (SDF) of various recurrence intervals for all stations are derived. Gumbel and BB1 copulas are typically chosen because they are more suitable for the DD and DS series, respectively. To be more precise, Joe copulas are only found over the basin at average elevation, but BB6 copulas can be seen at both high and low elevations. Furthermore, the BB7 copula is captured in low and high elevations. The results of the best-fitted copula are convenient with the Kiafar et al. (2020) and Nabaei et al. (2019) papers. The spatial distribution joint return period of T DS and T′ DS considering drought risk assessment based on percentile method (25%, 50%, 75%, and 95%) are displayed at four different categories (mild, moderate, severe, and extreme drought). Based on categorizing drought return periods, their interpretation provides clear information about drought analysis. Lower return periods should be considered high risk for drought occurrence that could negatively affect water quality, water suppliers, and soil moisture. Hence, when considering the SPI3 time scale, the high drought risk generally is expected to be experienced in the middle part of the basin. The lower drought risk is shown in the northern region for T DS , while a higher drought risk (lower return period) is detected under the mild and moderate drought categories for the T′ DS return period. Some researchers have studied on drought investigation on Ceyhan Basin (Gumus and Algin 2017;Altın et al. 2020;Dikiċi̇ and Aksel 2021;Topçu et al. Yuce et al. 2022). Yuce and Esit (2021) and Dikiċi̇ and Aksel (2021) observed to compare drought indices on the basin. However, these studies considered only one drought parameter (duration, severity, or intensity) instead of using multiple parameters. Unlike previous studies, a comprehensive drought analysis using copula function is investigated over the Ceyhan Basin, Turkey. In addition, results are also considered at different elevation levels. The finding of this paper may help the planning of future water resource management and mitigation of drought impact in the Ceyhan Basin.