Conditional Copula-Based Spatial–Temporal Drought Characteristics Analysis—A Case Study over Turkey

Abstract

In this study, commonly used copula functions belonging to Archimedean and Elliptical families are fitted to the univariate cumulative distribution functions (CDF) of the drought characteristics duration ($\text{LD}$), average severity ($\overline{\mathrm{S}}$), and average areal extent ($\overline{\mathrm{A}}$) of droughts obtained using standardized precipitation index (SPI) between 1960 and 2013 over Ankara, Turkey. Probabilistic modeling of drought characteristics with seven different fitted copula functions and their comparisons with independently estimated empirical joint distributions show normal copula links drought characteristics better than other copula functions. On average, droughts occur with an average $\text{LD}$ of 6.9 months, $\overline{\mathrm{S}}$ of 0.94, and $\overline{\mathrm{A}}$ of 73%, while such a drought event happens on average once in every 6.65 years. Results also show a very strong and statistically significant relation between $\overline{\mathrm{S}}$ and $\overline{\mathrm{A}}$, and drought return periods are more sensitive to the unconditioned drought characteristic, while return periods decrease by adding additional variables to the analysis (i.e., trivariate drought analysis compared to bivariate).

1. Introduction

Droughts are a climatic phenomenon that may impact large and small regions alike for long or short time periods and influence almost all aspects of society [1]. Historically, droughts have had negative effects on domestic and agricultural water supply needs [2] and consequently altered the behavior of ecosystems [3,4]. Plans, designs, and operational applications for water resources systems, therefore, require accurate estimation of different drought characteristics, like duration ($\text{LD}$), average severity ($\overline{\mathrm{S}}$), average areal extent ($\overline{\mathrm{A}}$), and re-occurrence rate under various scenarios.

Drought characteristics have been analyzed using various methods that have evolved over time. Earlier studies investigated the analytical or stochastic nature of past droughts using runs theory. Based on this theory, observed time series of drought-related variables are divided into wet and dry periods (positive and negative runs, respectively; see Figure 1) depending on a given threshold or critical level [5,6]. Later studies focused more on probabilistic drought estimation using ground station-based datasets [7,8,9]. The scarcity of ground stations and short recording lengths motivated the following studies to use synthetic datasets to analyze the frequencies of drought events [10,11]. Since the spatial–temporal relationships among drought characteristics are complex, recent studies have diversified their efforts by analyzing the return periods of droughts using various $\overline{\mathrm{S}}$ or $\overline{\mathrm{A}}$ scenarios [12,13,14].

The copula approach [15] provides an ideal test bed in which to analyze multivariate probabilistic problems without making explicit assumptions about the marginal or joint distributions of the variables involved in the estimation problem. Perhaps this ability has motivated studies to implement a multivariate probabilistic approach via copula functions to solve many different problems in hydrologic science and water resources-related applications [16,17,18,19].

Given that drought return periods can be immediately formulated in a probabilistic framework, copula-related methods have the necessary utility required to investigate different drought characteristics. Shiau [13] used bivariate copula functions to investigate bivariate drought analysis, including joint probabilities and return periods. Song and Singh [20] expanded this drought frequency analysis from bivariate to trivariate via a new method that constructs trivariate copulas to describe the joint distribution function of the temporal characteristics of meteorological drought events. Recently, Xu et al. [14] used a trivariate drought identification method within the copula approach to calculate the joint probability and return period of different drought spatial–temporal characteristic pairs. Overall, these studies demonstrated the advantages of the copula approach on bivariate and trivariate modeling of drought characteristics.

Certain characteristics of drought may precede others. For example, droughts with longer durations may tend to have higher maximum severity, or it may take longer to dissipate severe droughts than droughts with moderate severity. Such dependencies between drought characteristics can be immediately studied in a conditional copula framework, particularly when the occurrence of certain characteristics signals the occurrence of other characteristic functions in multivariate problems. Even though joint copula functions have been implemented in many studies, the implementation of conditional copulas remained limited in the literature [14,21]; hence there is still more room for more conditional copula-based applications to fully comprehend the utility of these functions.

Drought characteristics could be profoundly linked both in time and space. Accordingly, knowledge (i.e., observations) of some drought characteristics in time or space may help us predict others. This temporal and/or spatial dependency between different drought characteristics may immediately be exploited using copula functions, particularly using conditional probabilities. Even though spatial–temporal dependencies of drought characteristics have been investigated before using copula functions [22,23,24], such investigations have not been carried out in a conditional framework before.

The main objective of this study is to implement a copula-based methodology to analyze the joint and conditional dependency among various spatial–temporal characteristics of droughts, including $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$. The analysis utilizes univariate cumulative distribution functions (CDF) of drought characteristics obtained over Ankara, Turkey to establish joint and conditional relations between them, using bivariate and trivariate copula functions to calculate expected return periods of droughts with different characteristics scenarios.

2. Materials and Methods

In the following section, copula functions and their application to probabilistic drought characterizations are explained; later the use of bivariate and trivariate functions (i.e., the use of return periods) in drought analysis is elaborated.

2.1. Copula Functions

Copula functions link the univariate CDFs of random variables with their joint CDFs [13,15]. Estimations of these univariate CDFs are trivial and the sample datasets are often sufficient, while the analytical solution for the joint CDFs is not immediately clear. At this point, copula functions link univariate CDFs and create their multivariate (e.g., bivariate or trivariate) CDFs.

Based on Sklar’s theorem, bivariate and trivariate joint CDFs denoted as $\mathrm{F}(\mathrm{u},\mathrm{v})$ and $\mathrm{F}(\mathrm{u},\mathrm{v},\mathrm{w})$, respectively, can be obtained in the form:

$$\mathrm{F}(\mathrm{u},\mathrm{v})=\mathrm{C}(\mathrm{F}(\mathrm{u}),\mathrm{F}(\mathrm{v})),$$

(1)

$$\mathrm{F}(\mathrm{u},\mathrm{v},\mathrm{w})=\mathrm{C}(\mathrm{F}(\mathrm{u}),\mathrm{F}(\mathrm{v}),\mathrm{F}(\mathrm{w})),$$

(2)

where $\mathrm{u}$, $\mathrm{v}$, and $\mathrm{w}$ are random variables (e.g., $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$); $\mathrm{C}( )$ is the copula function; and univariate CDFs $\mathrm{F}(\mathrm{u})$, $\mathrm{F}(\mathrm{v})$, and $\mathrm{F}(\mathrm{w})$ are inputs to these function. Here $\mathrm{F}(\mathrm{u})$, $\mathrm{F}(\mathrm{v})$, and $\mathrm{F}(\mathrm{w})$ are the univariate CDFs of $\mathrm{u}$, $\mathrm{v}$, and $\mathrm{w}$ respectively, and they are defined as

$$\mathrm{F}(\mathrm{u})=\mathrm{P}(\mathrm{U}\le \mathrm{u}),$$

(3)

$$\mathrm{F}(\mathrm{v})=\mathrm{P}(\mathrm{V}\le \mathrm{v}),$$

(4)

$$\mathrm{F}(w)=\mathrm{P}(\mathrm{W}\le \mathrm{w}),$$

(5)

where $\mathrm{P}(\mathrm{U}\le \mathrm{u})$, $\mathrm{P}(\mathrm{V}\le \mathrm{v})$, and $\mathrm{P}(\mathrm{W}\le \mathrm{w})$ are the probabilities that random variables $\mathrm{U}$, $\mathrm{V}$, and $\mathrm{W}$ (i.e., different drought characteristics) takes values smaller than $\mathrm{u}$, $\mathrm{v}$, and w.

In Equations (1) and (2), the copula functions link the univariate CDF to bivariate and trivariate joint CDFs ($\mathrm{F}(\mathrm{u},\mathrm{v})$ and $\mathrm{F}(\mathrm{u},\mathrm{v},\mathrm{w})$). Using these joint CDFs, bivariate and trivariate conditional CDFs can be found by utilizing the multiplication rule:

$$\mathrm{F}(\mathrm{u}|\mathrm{v})=\frac{\mathrm{F}(\mathrm{u},\mathrm{v})}{\mathrm{F}(\mathrm{v})},$$

(6)

$$\mathrm{F}(\mathrm{u},\mathrm{v}|\mathrm{w})=\frac{\mathrm{F}(\mathrm{u},\mathrm{v},\mathrm{w})}{\mathrm{F}(\mathrm{w})}.$$

(7)

In Equation (7), $\mathrm{F}(\mathrm{u},\mathrm{v}|\mathrm{w})$ is given as an example of the trivariate CDF, while many other trivariate CDF (e.g., $\mathrm{F}(\mathrm{u}|\mathrm{v},\mathrm{w})$, $\mathrm{F}(\mathrm{v}|\mathrm{u},\mathrm{w})$, etc.) can be found in a similar fashion.

In this study, the bivariate (Equation (1)) and trivariate (Equation (2)) joint CDF of $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$ are obtained using seven widely used copula functions belonging to Archimedean and Elliptical copula families. These copula functions and their parameter spaces are shown in

Table 1. Equations of Copula functions. Here u and ν are two dependent univariate variables, df is the degree of freedom, $\mathsf{\theta }$ and $\mathsf{\rho }$ are Copula dependence parameters, $\varnothing$ is the CDF of standard univariate Gaussian distribution, and ${\mathrm{t}}_{\text{df}}$ is the t–student distribution function. AMH and GH refer to the Ali–Mikhail–Haq and Gumbel–Hougaard distributions, respectively.

Function (Family)	Joint Cumulative Distribution Function, $\mathbf{C}\mathbf{(}\mathbf{u}\mathbf{,}\mathbf{v}\mathbf{)}$	Parameter Range
AMH (Archimedean)	$\frac{\text{uv}}{1-\mathsf{\theta }(1-\mathrm{u})(1-\mathrm{v})}$	$-1\le \mathsf{\theta }\le 1$
Clayton (Archimedean)	${[\text{max}(0,{\mathrm{u}}^{-\mathsf{\theta }}+{\mathrm{v}}^{-\mathsf{\theta }}-1)]}^{-1/\mathsf{\theta }}$	$0\le \mathsf{\theta }$
Frank (Archimedean)	$\frac{-1}{\mathsf{\theta }}\text{ln}(1+\frac{({\mathrm{e}}^{-\mathsf{\theta }\mathrm{u}}-1)({\mathrm{e}}^{-\mathsf{\theta }\mathrm{v}}-1)}{{\mathrm{e}}^{-\mathsf{\theta }}-1})$	$\mathsf{\theta }\ne 0$
GH (Archimedean)	${\mathrm{e}}^{{[{(-\ln \mathrm{u})}^{\mathsf{\theta }}+{(-\ln \mathrm{v})}^{\mathsf{\theta }}]}^{\frac{1}{\mathsf{\theta }}}}$	$1\le \mathsf{\theta }$
Joe (Archimedean)	$1-{({(1-\mathrm{u})}^{\mathsf{\theta }}+{(1-\mathrm{v})}^{\mathsf{\theta }}-{(1-\mathrm{u})}^{\mathsf{\theta }}{(1-\mathrm{v})}^{\mathsf{\theta }})}^{\frac{1}{\mathsf{\theta }}}$	$1\le \mathsf{\theta }$
Normal (Elliptical)	${\displaystyle \underset{-\infty }{\overset{{\varnothing }^{-1}(\mathrm{u})}{\int }}\underset{-\infty }{\overset{{\varnothing }^{-1}(\mathrm{v})}{\int }}}\frac{1}{2\mathsf{\pi }{(1-{\mathsf{\rho }}^2)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\mathrm{e}}^{-\frac{{\mathrm{u}}^2-2\mathsf{\rho }\text{uv}+{\mathrm{v}}^2}{2(1-{\mathsf{\rho }}^2)}}\text{dudv}$	$-1\le \mathsf{\rho }\le 1$
t–student (Elliptical)	${\displaystyle \underset{-\infty }{\overset{{{\mathrm{t}}_{\text{df}}}^{-1}(\mathrm{u})}{\int }}\underset{-\infty }{\overset{{{\mathrm{t}}_{\text{df}}}^{-1}(\mathrm{v})}{\int }}}\frac{1}{2\mathsf{\pi }{(1-{\mathsf{\rho }}^2)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{(1+\frac{{\mathrm{u}}^2-2\mathsf{\rho }\text{uv}+{\mathrm{v}}^2}{\text{df}(1-{\mathsf{\rho }}^2)})}^{-\frac{\text{df}+2}{2}}\text{dudv}$	$-1\le \mathsf{\rho }\le 1$ $\text{df}\ge 1$

. For further details about these copula functions, see the studies of Joe [25], Shiau [13], and Nelsen [26].

The required copula dependence parameters ($\mathsf{\theta }$ and $\mathsf{\rho }$) for the seven copula functions are separately estimated for each copula function using the inference function for margins (IFM) method [25]. This method involves estimating $\mathrm{F}(\mathrm{u})$ and $\mathrm{F}(\mathrm{v})$ for the variables and then estimating the $\mathsf{\theta }$ or $\mathsf{\rho }$ parameter that maximizes the Kendall’s tau correlation using these univariate CDFs. These obtained parameters are later used to make copula function-based joint CDF predictions, which are referred to as “theoretical joint CDF estimation” ($\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{t}}$) in this study.

The above described copula function-based $\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{t}}$ values are later independently validated using “empirical joint CDF estimates” ($\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{e}}$) obtained using rank-based joint distribution estimates obtained via the copula package [27] in the R programming language [28]. Datasets are split for parameter estimation and validation with the leave-one-out cross-validation method to avoid over-fitting. Validation efforts include calculation of root mean square error (RMSE) using $\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{t}}$ and $\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{e}}$:

$$\text{RMSE}=\sqrt{\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{t}}-\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{e}})}^2}{\mathrm{n}}},$$

(8)

where n is the length of the drought characteristic datasets. Among the seven copula functions, the one that yields the lowest RMSE value is later selected for further estimation of bivariate and trivariate return periods.

2.2. Return Periods

The expected value of the drought inter-arrival time ($\text{DIT}$; the time between two successive drought events) is one of the most commonly estimated drought characteristics in hydrological and water resources systems applications. $\mathrm{T}$ may be estimated for univariate (e.g., $\text{LD}\ge 10$ months) as well as multivariate (e.g., $\text{LD}\ge 10$ & $\overline{\mathrm{S}}\ge 1.0$; or $\text{LD}\ge 10$|$\overline{\mathrm{S}}\ge 1.0$) scenarios. Shiau and Shen [29] theoretically derived $\mathrm{T}$ for drought events with univariate drought characteristics (e.g., $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$) as a function of the expected $\text{DIT}$ in the form

$${\mathrm{T}}_{\mathrm{U}\ge \mathrm{u}}=\frac{\mathrm{E}(\text{DIT})}{\mathrm{P}(\mathrm{U}\ge \mathrm{u})},$$

(9)

where $\mathrm{E}(\text{DIT})$ is the expected value of $\text{DIT}$; and ${\mathrm{T}}_{\mathrm{U}\ge \mathrm{u}}$ is the return period of a drought event, defined by the variable $\mathrm{u}$. $\mathrm{E}(\text{DIT})$ is calculated as the ratio of the duration of entire study to the total number of droughts. Following Shiau [13], it is also possible to expand this univariate $\mathrm{T}$ estimation analysis to a bivariate case by characterizing events with bivariate joint or conditional probabilities in the forms:

$${\mathrm{T}}_{(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v})}=\frac{\mathrm{E}(\text{DIT})}{\mathrm{P}(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v})}$$

(10)

$${\mathrm{T}}_{(\mathrm{U}\ge \mathrm{u} \text{or} \mathrm{V}\ge \mathrm{v})}=\frac{\mathrm{E}(\text{DIT})}{\mathrm{P}(\mathrm{U}\ge \mathrm{u} \text{or} \mathrm{V}\ge \mathrm{v})}$$

(11)

$${\mathrm{T}}_{(\mathrm{U}\ge \mathrm{u}|\mathrm{V}\ge \mathrm{v})}=\frac{\mathrm{E}(\text{DIT})}{\mathrm{P}(\mathrm{U}\ge \mathrm{u},|\mathrm{V}\ge \mathrm{v})}.$$

(12)

Similarly, the trivariate case can be obtained in the forms:

$${\mathrm{T}}_{(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w})}=\frac{\mathrm{E}(\text{DIT})}{\mathrm{P}(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w})}$$

(13)

$${\mathrm{T}}_{(\mathrm{U}\ge \mathrm{u} \text{or} \mathrm{V}\ge \mathrm{v} \text{or} \mathrm{W}\ge \mathrm{w})}=\frac{\mathrm{E}(\text{DIT})}{\mathrm{P}(\mathrm{U}\ge \mathrm{u} \text{or} \mathrm{V}\ge \mathrm{v} \text{or} \mathrm{W}\ge \mathrm{w})}$$

(14)

$${\mathrm{T}}_{(\mathrm{U}\ge \mathrm{u}|(\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w}))}=\frac{\mathrm{E}(\text{DIT})}{\mathrm{P}(\mathrm{U}\ge \mathrm{u}|(\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w}))},$$

(15)

where the $\mathrm{P}(\dots )$ values in the denominators of Equations (10)–(12) are the bivariate, and of Equations (13)–(15) are the multivariate joint probabilities of drought characteristics $\mathrm{u}$, $\mathrm{v}$, and $\mathrm{w}$; and ${\mathrm{T}}_{(\dots )}$ values on the left-hand side of Equations (10)–(15) are the corresponding return periods. By definition, the $\mathrm{P}(\dots )$ terms in Equations (9)–(15) can be found using the CDFs of the related drought characteristics. Below are some examples of univariate, bivariate, and trivariate cases:

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u})=1-{\mathrm{F}}_{\mathrm{U}}(\mathrm{u})$$

(16)

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u}|\mathrm{V}\ge \mathrm{v})=1-{\mathrm{F}}_{\mathrm{U},\mathrm{V}}(\mathrm{u}|\mathrm{v})$$

(17)

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w})=1-{\mathrm{F}}_{\mathrm{U},\mathrm{V},\mathrm{W}}(\mathrm{u},\mathrm{v},\mathrm{w}).$$

(18)

Accordingly, the bivariate probability $\mathrm{P}(\dots )$ terms in Equations (10)–(12) can be found using the univariate CDF estimates and the bivariate copula functions in the forms:

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v})=1-{\mathrm{F}}_{\mathrm{U}}(\mathrm{u})-{\mathrm{F}}_{\mathrm{V}}(\mathrm{v})+\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{V}}(\mathrm{v}))$$

(19)

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u} \text{or} \mathrm{V}\ge \mathrm{v})=1-\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{V}}(\mathrm{v}))$$

(20)

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u}|\mathrm{V}\ge \mathrm{v})=\frac{\mathrm{P}(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v})}{\mathrm{P}(\mathrm{V}\ge \mathrm{v})}=\frac{1-{\mathrm{F}}_{\mathrm{U}}(\mathrm{u})-{\mathrm{F}}_{\mathrm{V}}(\mathrm{v})+\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{V}}(\mathrm{v}))}{1-{\mathrm{F}}_{\mathrm{V}}(\mathrm{v})}.$$

(21)

The trivariate probabilities in Equations (13)–(15) can be found using the univariate CDF estimates, and bivariate and trivariate copula functions. Below is an example of a trivariate conditional case, while many other forms could be obtained in a similar fashion:

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u}|(\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w}))=\frac{\mathrm{P}(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w})}{\mathrm{P}(\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w})},$$

(22)

where the bivariate and trivariate joint probability terms on the right-hand side can be expressed in terms of copula functions. Below are some examples of such expressions linking multivariate probabilistic terms with copula functions:

$$\mathrm{P}(\mathrm{V}\le \mathrm{v},\mathrm{W}\le \mathrm{w})=\mathrm{C}({\mathrm{F}}_{\mathrm{V}}(\mathrm{v}),{\mathrm{F}}_{\mathrm{W}}(\mathrm{w}))$$

(23)

$$\mathrm{P}(\mathrm{V}\ge \mathrm{v},\mathrm{W}\ge \mathrm{w})=1-{\mathrm{F}}_{\mathrm{V}}(\mathrm{v})-{\mathrm{F}}_{\mathrm{W}}(\mathrm{w})+\mathrm{C}({\mathrm{F}}_{\mathrm{V}}(\mathrm{v}),{\mathrm{F}}_{\mathrm{W}}(\mathrm{w}))$$

(24)

$$\mathrm{P}(\mathrm{U}\le \mathrm{u},\mathrm{V}\le \mathrm{v},\mathrm{W}\le \mathrm{w})=\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{V}}(\mathrm{v}),{\mathrm{F}}_{\mathrm{W}}(\mathrm{w}))$$

(25)

$$\begin{array}{cc}\hfill \mathrm{P}(\mathrm{U}\ge \mathrm{u},\mathrm{V}\ge \mathrm{v},& \mathrm{W}\ge \mathrm{w})\hfill \\ & =1-{\mathrm{F}}_{\mathrm{U}}(\mathrm{u})-{\mathrm{F}}_{\mathrm{V}}(\mathrm{v})-{\mathrm{F}}_{\mathrm{W}}(\mathrm{w})+\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{V}}(\mathrm{v}))+\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{W}}(\mathrm{w}))\hfill \\ & +\mathrm{C}({\mathrm{F}}_{\mathrm{V}}(\mathrm{v}),{\mathrm{F}}_{\mathrm{W}}(\mathrm{w}))-\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{V}}(\mathrm{v}),{\mathrm{F}}_{\mathrm{W}}(\mathrm{w}))\hfill \end{array}$$

(26)

$$\mathrm{P}(\mathrm{U}\ge \mathrm{u} \text{or} \mathrm{V}\ge \mathrm{v} \text{or} \mathrm{W}\ge \mathrm{w})=1-\mathrm{C}({\mathrm{F}}_{\mathrm{U}}(\mathrm{u}),{\mathrm{F}}_{\mathrm{V}}(\mathrm{v}),{\mathrm{F}}_{\mathrm{W}}(\mathrm{w})).$$

(27)

Equations (9)–(15) describe how $\mathrm{T}$ values are estimated for univariate (9), bivariate (10)–(12), and trivariate (13)–(15) cases. These equations are flexible enough to accommodate investigations of various spatial–temporal drought characteristics like $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$ (e.g., $\mathrm{T}$ of a drought event with $\text{LD}\ge 10|\overline{\mathrm{S}}>0.80 \& \overline{\mathrm{A}}>0.85$ can be calculated). Once the multivariate joint probabilities of the drought characteristics are determined through these copula functions, $\mathrm{T}$ for these drought characteristics can be determined by plugging these probability estimates into Equations (10)–(15).

2.3. Standardized Precipitation Index (SPI)

In this study, SPI [30] is selected to analyze drought events. It is worth noting that this selection is arbitrary and inconsequential, while other drought indices (e.g., Percent of Normal, Palmer Drought Severity Index, Crop Moisture Index, Surface Water Supply Index, Reclamation Drought Index, etc.) could have been used as well.

SPI can be calculated for different time scales (i.e., 3, 6, 12, etc.), while the obtained time series (SPI-3, SPI-6, SPI-12, etc.) are used as drought indicators in applications with different goals. From a drought management point of view, managers may prefer SPI-3 as an indicator to monitor short-term anomalies, SPI-12 to monitor persistent events, or SPI-6 to have a flavor of both short-term anomalies and persistent events [31]. Accordingly, in this study SPI-6 is selected for the drought analysis to quantify the deficits of precipitation. For this purpose, all of the monthly precipitation time series are stored in a matrix ${\mathrm{P}}_{\mathrm{s},\mathrm{m}}$, where $\mathrm{s}$ refers to each station and $\mathrm{m}$ refers to each month. Later six-month moving-window precipitation sums (${\text{PT}}_{\mathrm{s},\mathrm{m}}$) are obtained in the form:

$${\text{PT}}_{\mathrm{s},\mathrm{m}}={\Sigma }_{m-5}^m{\mathrm{P}}_{\mathrm{s},\mathrm{m}}.$$

(28)

These ${\text{PT}}_{\mathrm{s},\mathrm{m}}$ values are spatially averaged to obtain representative values for the entire study area in the form:

$${\text{PT}}_{\mathrm{m}}={\sum }_{s=1}^D{\mathrm{P}}_{\mathrm{s},\mathrm{m}},$$

(29)

where $\mathrm{D}$ is the total number of stations. Following McKee et al. [30], these ${\text{PT}}_{\mathrm{s},\mathrm{m}}$ and ${\text{PT}}_{\mathrm{m}}$ time series are later fitted to gamma distribution, where the probability of any given precipitation value is obtained using this fitted gamma distribution. These probability estimates obtained using ${\text{PT}}_{\mathrm{s},\mathrm{m}}$ and ${\text{PT}}_{\mathrm{m}}$ are later converted into ${\text{SPI}}_{\mathrm{s},\mathrm{m}}$ and ${\text{SPI}}_{\mathrm{m}}$ values, respectively, using the inverse of the standard normal distribution. Here, ${\text{SPI}}_{\mathrm{s},\mathrm{m}}$ values are calculated for each station and month separately, while ${\text{SPI}}_{\mathrm{m}}$ is representative for the entire study area and calculated for each month.

2.4. Drought Characteristics

$\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$ are the three drought characteristics considered in this study. The definition of the dry and the non-dry periods during the study period is made using ${\text{SPI}}_{\mathrm{m}}$ time series. Following McKee et al. [30], the dry periods are defined to have continuously negative SPI values and have at least one SPI value below −1.0 (e.g., the time series (−0.23, −0.76, −1.15, −0.52, −0.05) constitute a drought event, while the time series (−0.23, −0.76, −0.99, −0.52, −0.05) does not). Accordingly, the periods that are not considered as dry are considered “non-dry” periods. $\text{LD}$ and the length of non-dry periods ($\text{LnD}$) values are calculated as:

$$\text{LD}=({\mathrm{t}}_{\text{de}}-{\mathrm{t}}_{\text{ds}})+1$$

(30)

$$\text{LnD}=({\mathrm{t}}_{\text{nde}}-{\mathrm{t}}_{\text{nds}})+1,$$

(31)

where ${\mathrm{t}}_{\text{ds}}$ and ${\mathrm{t}}_{\text{de}}$ are the months a particular drought period started and ended, respectively; and ${\mathrm{t}}_{\text{nds}}$ and ${\mathrm{t}}_{\text{nde}}$ are the months a particular non-dry period started and ended, respectively. $\text{DIT}$ is obtained as

$$\text{DIT}=\text{LD}+\text{LnD}.$$

(32)

Even though ${\text{SPI}}_{\mathrm{s},\mathrm{m}}$ and ${\text{SPI}}_{\mathrm{m}}$ values are available for each month, drought characteristic values ($\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$) are calculated for each drought period. Among these drought characteristics, $\text{LD}$ values are calculated using ${\text{SPI}}_{\mathrm{m}}$ time series. For each drought period, $\overline{\mathrm{S}}$ is calculated as:

$$\overline{\mathrm{S}}=|\frac{1}{\text{LD}}{\sum }_{m={\mathrm{t}}_{\text{ds}}}^{t_{\text{de}}}{\text{SPI}}_{\mathrm{m}}|.$$

(33)

Separate drought definitions are made for the region and the stations using ${\text{SPI}}_{\mathrm{m}}$ and ${\text{SPI}}_{\mathrm{s},\mathrm{m}}$ values, respectively. If the given region is under drought conditions, then individual stations are assessed as to whether or not they individually satisfy the drought conditions using ${\text{SPI}}_{\mathrm{s},\mathrm{m}}$ values. Then the number of months of each station that are under drought condition during the regional drought period is calculated. The fraction of station-based drought conditions compared to the duration of regional drought length is defined as $\overline{\mathrm{A}}$. If the $\text{SPI}$ values for a particular month are only available for four stations, then $\overline{\mathrm{A}}$ is calculated as the fraction of these four stations that are under drought considering the length of the particular drought event. Once drought characteristics $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$ are obtained for each drought period separately, the univariate, bivariate, and trivariate $\mathrm{T}$ values of these characteristics are calculated using Equations (10)–(15).

2.5. Drought Frequency and Visual Analysis

The visual presentation of the multivariate CDFs of the drought characteristics is instructive for building knowledge about the drought events. However, such analysis requires long datasets for visually presentable plot generation. In the absence of a sufficiently long time series, it is trivial to synthetically generate sufficient datasets, while such simulations require the distribution of the particular drought characteristic ($\text{LD}$, $\overline{\mathrm{S}}$, or $\overline{\mathrm{A}}$) to be known in advance. In many cases, observed data do not follow a specific distribution, and their distributions likely vary depending on the variables and time series [21]. In this study, $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$ are respectively fitted to various distributions (the seven advanced with variable support, 23 non-negative, and eight bounded available distributions;

Table 2. Fitted distribution functions (total of 38 types).

Advanced with Variable Support	Non-Negative	Non-Negative	Non-Negative	Bounded
Gen. extr. val.	Burr	Gamma	Pareto	Beta
Gen. logistic	Chi-squared	Gen. gamma	Pareto 2	Johnson SB
Gen. Pareto	Dagum	Inv. Gaussian	Pearson 5	Kumaraswamy
Log–Pearson 3	Erlang	Levy	Pearson 6	Pert
Phased Bi–exp.	Exp.	Log-gamma	Rayleigh	Power func.
Phased bi–Weibull	F	Log-logistic	Rice	Reciprocal
Wakeby	Fatigue life	Lognormal	Weibull	Triangular
	Frechet	Nakagami		Uniform

) in the EasyFit program [32] to find the distribution that resulted in the best fit to the CDF of the observed drought characteristics data. The best CDF fit is separately found for each drought characteristic using the Maximum Likelihood method [33] by optimizing the Kolmogrov–Smirnov index [34] within the EasyFit program. The distributions that resulted in the best CDF fits are later used to synthetically generate sufficiently long datasets. Here it is stressed that these synthetically generated datasets are only used for visual representations but not in other analysis.

2.6. Study Area and the Datasets

Analysis of drought events is very critical over arid and semi-arid regions, like Ankara, Turkey, where water management projects require knowledge of the water status during a particular drought period. Among the different types of droughts (meteorological, hydrological, agricultural, and socioeconomic), this study focuses on the analysis of the spatial–temporal characteristics of meteorological droughts using monthly precipitation records obtained over six meteorological stations (Beypazari, Esenboga, Kecioren, Kizilcahamam, Nallihan, and Polatli) distributed over Ankara, Turkey (Figure 2). These stations are maintained and the datasets are distributed by The General Directorate of the Turkish State Meteorological Service (TSMS). Precipitation datasets over these six stations between January 1960 and February 2014 are obtained from TSMS. Precipitation time series for the first 60 months and the last two months for Nallihan and Esenboga stations, respectively (out of a total of 650 months). These missing data are ignored for these missing periods, while they constitute only <2% of the total datasets. More information about these stations can be found in the study of Sonmez [35].

3. Results and Discussion

The drought characteristics $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$ were calculated for each drought period separately (

Table 3. Drought characteristics defied by using SPI series values.

Drought	Start of Drought		End of Drought		LD	$\overline{\mathbf{S}}$	$\overline{\mathbf{A}}$
	Year	Month	Year	Month
1	1960	9	1961	1	5	0.9	0.8
2	1961	8	1962	1	6	0.68	0.67
3	1962	7	1962	11	5	0.94	0.76
4	1965	8	1965	12	5	1.16	0.92
5	1966	9	1967	2	6	0.82	0.78
6	1967	9	1968	1	5	1.05	0.87
7	1969	9	1969	12	4	1.28	0.96
8	1970	8	1971	2	7	0.75	0.64
9	1973	3	1974	4	14	0.84	0.82
10	1974	11	1975	1	3	0.84	0.44
11	1976	7	1976	12	6	0.39	0.17
12	1977	6	1978	1	8	1.03	0.81
13	1978	8	1978	12	5	0.79	0.67
14	1979	7	1979	12	6	1.07	0.5
15	1980	8	1980	12	5	0.95	0.9
16	1981	8	1981	11	4	0.79	0.5
17	1982	10	1983	5	8	0.93	0.94
18	1984	9	1985	3	7	1.28	1
19	1985	7	1986	1	7	1.05	0.74
20	1986	7	1986	12	6	0.86	0.81
21	1987	7	1988	2	8	0.91	0.88
22	1988	9	1989	10	14	0.82	0.92
23	1991	8	1992	2	7	0.53	0.48
24	1992	6	1993	2	9	0.52	0.48
25	1993	7	1995	2	20	1.05	0.86
26	1996	8	1996	12	5	0.48	0.33
27	1998	10	1999	1	4	0.68	0.33
28	2000	9	2001	4	8	1.28	0.94
29	2001	8	2001	11	4	0.76	0.67
30	2003	7	2003	12	6	1.43	0.97
31	2004	7	2005	2	8	1.34	0.98
32	2005	10	2006	1	4	0.78	0.46
33	2007	3	2008	2	12	1.08	0.86
34	2008	6	2008	12	7	0.86	0.48
35	2009	9	2009	12	4	0.93	0.63
36	2012	7	2012	12	6	1.33	1
37	2013	7	2014	2	8	1.46	0.91
Statistics of the above given 37 drought periods				Min	3	0.39	0.17
				Mean	6.92	0.94	0.73
				Max	20	1.46	1

). Accordingly, higher $\overline{\mathrm{S}}$ values indicate more severe drought conditions and values closer to 0 indicate near to normal conditions. For example, the drought event starting in July 1985 and ending in January 1986 (Table 3) has an $\text{LD}$ of seven months, $\overline{\mathrm{S}}$ of 1.05, and $\overline{\mathrm{A}}$ of 0.74. Here $\overline{\mathrm{S}}$ is calculated as the absolute value of the average of seven SPI values: $|-(0.28+1.04+1.50+1.42+1.60+1.33+0.17)/7|=1.05$; and $\overline{\mathrm{A}}$ is calculated as the fraction of total months under drought condition (i.e., the regional average SPI time series show the drought duration is seven months, while the six separate SPI time series over each station show 31 months are under drought conditions out of a total of six stations × seven months = 42 months, implying $\overline{\mathrm{A}}=31/42=0.74$). Overall results show that drought conditions influenced the study area during 39% of the total study period (summation of all $\text{LD}$/650 months), while these drought events, on average, have $\overline{\mathrm{S}}$ of 0.94, last for 6.92 months, and impact 73% of the area covered by the stations (Table 3).

Droughts tend to start during drier months (July–September) and end during wetter winter months (November–February). This result is expected as SPI-6 considers six-monthly accumulated precipitation values (Equation (28)): SPI-6 values for July–September are calculated using relatively drier months and values for November–February are calculated using relatively wetter months, where the summer and the winter months are climatologically dry and wet, respectively. This systematic pattern is consistent with the drought definition of many water resources managers who are interested in the lack of water. Removal of this strong SPI-6 seasonality may result in reduced SPI-6 sensitivity to actual water deficit (i.e., anomalies obtained during drier months may not necessarily mean the same “water deficit” during wetter months), hence SPI-6 seasonality is retained in this study.

Strong linear dependence is found between $\overline{\mathrm{S}}$ and $\overline{\mathrm{A}}$ (correlations are statistically significant at 95% confidence level), while the relationships between other drought characteristics are much weaker (

Table 4. Kendall’s tau and Pearson’s correlations between drought characteristics LD, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$.

	$\mathbf{\text{LD}}\mathbf{–}\overline{\mathbf{S}}$	$\mathbf{\text{LD}}\mathbf{–}\overline{\mathbf{A}}$	$\overline{\mathbf{S}}\mathbf{–}\overline{\mathbf{A}}$
Kendall’s tau	0.15	0.27	0.63
Pearson‘s	0.10	0.32	0.79

). This implies that predictability studies involving $\overline{\mathrm{S}}$ and $\overline{\mathrm{A}}$ may contain high predictive skill, while particularly $\text{LD}$ prediction using $\text{LD}–\overline{\mathrm{S}}$ or $\text{LD}–\overline{\mathrm{A}}$ relations may not yield accurate results. Here, drought characteristics $\text{LD}$ and $\overline{\mathrm{A}}$ have boundaries (0 and 1.0, respectively); however, there are not many values close to these boundary values, hence the use of correlation statistics may not be problematic assuming these drought characteristics are approximately normal.

The univariate CDFs of drought characteristics are necessary to fit the copula functions [$\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{t}}$], so the best fitting copula function can be used for further analysis. Here, the best-fitting copula is obtained after performing validation efforts utilizing $\mathrm{C}{(\mathrm{u},\mathrm{v})}_{\mathrm{e}}$ as the truth. For the purpose of fitting copulas, univariate empirical CDF values are used while the theoretical univariate CDF values are used in the synthetic simulations described below (Figure 3). Among copula functions, by definition, the nature of the AMH requires Kendall’s tau between variables to be bounded between −0.18 and 0.33 [27]. Since the $\overline{\mathrm{S}}$–$\overline{\mathrm{A}}$ pair does not satisfy this requirement (Table 4), an AMH function could not be fitted to this pair. Leave-one-out type validation of copula parameters show that the elliptical copula family results in better performance in bivariate and trivariate joint CDF estimation analysis than other copula families, while normal copula seems to have the best overall performance (

Table 5. Parameters ($\mathsf{\theta }$ or $\mathsf{\rho }$) and RMSE values of bivariate and trivariate copula fitting simulations, where the three smallest RMSE values are written in bold for each analysis.

Copula Function	Variables of Bivariate Copulas						Variables of Trivariate Copula
	LD & $\overline{\mathbf{S}}$		LD & $\overline{\mathbf{A}}$		$\overline{\mathbf{S}}$ & $\overline{\mathbf{A}}$		LD & $\overline{\mathbf{S}}$ & $\overline{\mathbf{A}}$
	($\mathsf{\theta }$/$\mathsf{\rho }$)	RMSE	($\mathsf{\theta }$/$\mathsf{\rho }$)	RMSE	($\mathsf{\theta }$/$\mathsf{\rho }$)	RMSE	($\mathsf{\theta }$/$\mathsf{\rho }$)	RMSE
AMH	0.514	0.062	0.895	0.056	-	-	-	-
Clayton	0.280	0.064	0.728	0.058	2.395	0.043	0.930	0.067
Frank	1.272	0.061	2.425	0.058	8.170	0.030	3.026	0.065
GH	1.115	0.065	1.271	0.064	2.514	0.029	1.425	0.070
Joe	1.120	0.070	1.269	0.075	2.966	0.042	1.538	0.084
Normal	0.242	0.059	0.426	0.055	0.839	0.026	0.528	0.054
t-Student	0.244	0.061	0.426	0.056	0.838	0.026	0.528	0.056

).

Analysis of drought period counts show a drought event is expected on average every 18 months (i.e., $\mathrm{E}(\text{DIT})=650 \text{months}/37 \text{droughts}$). Univariate drought $\mathrm{T}$ calculations using the univariate CDF of variables (Figure 3) show on average drought events with $\text{LD}\ge 6.92$ are expected every 4.15 years, $\overline{\mathrm{S}}\ge 0.94$ every 3.06 years, and $\overline{\mathrm{A}}\ge 0.73$ every 2.53 years (

Table 6. Univariate drought return periods (T).

T (Years)	Drought Characteristics
	LD (Month)	$\overline{\mathbf{S}}$ (-)	$\overline{\mathbf{A}}$ (-)
1.46	3.00	-	-
1.47	3.24	-	0.17
1.48	3.26	0.39	0.18
2.53	5.71	0.87	0.73
3.06	6.21	0.94	0.79
4.15	6.92	1.02	0.86
40.57	14.20	1.44	1.00
48.19	15.01	1.46	-
116.83	20.00	-	-

). However, in many cases these univariate $\mathrm{T}$ estimates may not be sufficient for various drought investigations requiring bivariate or trivariate $\mathrm{T}$ information. These multivariate drought $\mathrm{T}$ values are calculated using the best fitting copula functions for the bivariate and the trivariate cases. As an example, $\mathrm{T}$ is calculated for the bivariate case as 7.13 years for the drought event with $\overline{\mathrm{S}}\ge 0.94$ and $\text{LD}\ge 6.92$; as 3.41 years for a drought event with $\overline{\mathrm{A}}\ge 0.73$ and $\overline{\mathrm{S}}\ge 0.94$; and for the trivariate case $\mathrm{T}$ is calculated as 6.65 years for a drought event with $\overline{\mathrm{S}}\ge 0.94$, $\overline{\mathrm{A}}\ge 0.73$, and $\text{LD}\ge 6.92$.

The examples given above could be reproduced for infinitely many drought characteristic pairs or triplets. However, plots showing their multivariate variations often are necessary to understand the exact nature of their multivariate relations. However, to generate such plots the lengths of the available 37 drought events are not sufficient. Instead, the drought characteristics are fitted to different distributions and the best fitting distributions are later used to generate synthetic drought characteristic datasets that are used to generate multivariate $\mathrm{T}$ plots (Figure 4 and Figure 5). Parameters result in the smallest Kolmogorov–Smirnov statistics (

Table 7. Kolmogorov–Smirnov statistics and parameters of best theoretical CDF distributions (Figure 3) for the drought characteristics $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$.

Drought Characteristic	Distribution	Parameters	Test Statistic
LD	Wakeby	α = 6.575; β = 3.344; γ = 1.470; δ = 0.34; ξ = 3.176	0.10
$\overline{\mathrm{S}}$	Pearson 5	α = 167.61; β = 563.17; γ = −2.44	0.08
$\overline{\mathrm{A}}$	Johnson SB	γ = −0.733; δ = 0.711; λ = 0.937; ξ = 0.088	0.08

) are used to generate the relevant synthetic datasets. Results show log-logistic, normal, and logistic distributions result in the best fit (i.e., lowest errors) for the drought characteristics of $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$, respectively (Table 7). Hence, these distributions are used to create synthetic $\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$ data.

Three drought characteristics ($\text{LD}$, $\overline{\mathrm{S}}$, and $\overline{\mathrm{A}}$) are considered in this study, hence three bivariate drought characteristic pairs ($\overline{\mathrm{S}}–\text{LD}$, $\overline{\mathrm{A}}–\text{LD}$, and $\overline{\mathrm{A}}–\overline{\mathrm{S}}$) are analyzed. Additionally, three bivariate drought operators (“and”, “or”, “conditional”) that are suited for various drought applications are also considered. Accordingly, a total of nine bivariate $\mathrm{T}$ maps were obtained for visual presentation efforts (Figure 4). In Figure 4, the left, middle, and right columns show “and”, “or”, and “conditional” scenarios, respectively; the top, middle, and bottom rows show the $\overline{\mathrm{S}}–\text{LD}$, $\overline{\mathrm{A}}–\text{LD}$, and $\overline{\mathrm{A}}$–S pairs, respectively.

The “and” scenario shows that the bivariate $\mathrm{T}$ has a higher sensitivity to $\text{LD}$ than $\overline{\mathrm{S}}$ or $\overline{\mathrm{A}}$: When $\text{LD}$ is included in the bivariate $\mathrm{T}$ maps (left middle and top panels), $\mathrm{T}$ values are more impacted by $\text{LD}$ than $\overline{\mathrm{A}}$ or $\overline{\mathrm{S}}$, particularly for higher values of $\text{LD}$, while the sensitivity difference between $\overline{\mathrm{A}}–\overline{\mathrm{S}}$ pairs is marginal. If a line is drawn from the lower left corner to the upper right corner of the $\overline{\mathrm{A}}–\overline{\mathrm{S}}$ bivariate $\mathrm{T}$ figure (lower-left panel of Figure 4), then the area to the left of this line shows the region where $\mathrm{T}$ is more sensitive to changes in $\overline{\mathrm{A}}$ than $\overline{\mathrm{S}}$ while the area on the right side of the line shows the region where $\mathrm{T}$ is more sensitive to changes in $\overline{\mathrm{S}}$ than $\overline{\mathrm{A}}$. These results suggest that the overall sensitivity of $\mathrm{T}$ to these drought characteristics can be compared via $\text{LD}$ > $\overline{\mathrm{A}}$–$\overline{\mathrm{S}}$. On average, bivariate $\mathrm{T}$ values are much higher, particularly for higher values of $\text{LD}$ ($\text{LD}>10$). Similar to the “and” case, the “or” case also show pareto front type plots but the pareto front direction is reversed (i.e., the middle column): on average almost one in 2.33 years a mild drought event with “$\overline{\mathrm{A}}\ge 0.73$ or $\overline{\mathrm{S}}\ge 0.94$” or “$\text{LD}\ge 6.91 \text{or} \overline{\mathrm{S}}\ge 0.94$” is observed. As expected, the average $\mathrm{T}$ value for the “or” case is much lower than for the “and” case as the probability of an “or” case is higher than an “and” case. In the case of “conditional” drought $\mathrm{T}$ plots, the $\mathrm{T}$ values are most sensitive to the unconditioned value rather than the given observation (i.e., the lines are vertical), while this sensitivity decreases as the drought impact increases (higher $\text{LD}$, $\overline{\mathrm{S}}$, or $\overline{\mathrm{A}}$). On average, the dynamic range of $\mathrm{T}$ values for conditional case are more similar to the “or” case, particularly for less severe drought events (e.g., for the bivariate case they range between 1.5 and 42 years). These conditional $\mathrm{T}$ values could be particularly useful if certain observations could be obtained. For example, if $\overline{\mathrm{S}}\ge 0.94$, then $\mathrm{T}$ values range between 1.86 and 9.95 years for $5\le \text{LD}\le 10$ (Figure 4, top right panel).

Figure 5 shows that the trivariate $\mathrm{T}$ values are conditioned on both $\overline{\mathrm{S}}$ and $\overline{\mathrm{A}}$. Given that it would be harder to digest 3D plots, trivariate cases are plotted for four separate $\overline{\mathrm{A}}$ values (0.65, 0.75, 0.85, and 0.95). Similar to conditional bivariate plots, conditional trivariate plots also show that $\mathrm{T}$ values are sensitive to the unconditioned variables (i.e., the plots are perpendicular to the unconditioned variable). Additionally, the conditional trivariate contour plots show that when additional observations are available compared to the bivariate case, then $\mathrm{T}$ values decrease (i.e., if $\overline{\mathrm{A}}\ge 0.65$, it is more likely that drought events with longer duration will happen). For example, the bivariate conditional plot for $\text{LD}|\overline{\mathrm{S}}\ge 0.7$ shows $\mathrm{T}$ values ranging between 1.95 and 42.04, while the trivariate conditional plot for $\text{LD}|\overline{\mathrm{S}}\ge 0.7,\overline{\mathrm{A}}\ge 0.95$ shows $\mathrm{T}$ values ranging between 2.05 and 13.69. The conditional $\mathrm{T}$ values have much higher sensitivity to $\text{LD}$ when $\overline{\mathrm{A}}\ge 0.75$ than when $\overline{\mathrm{A}}\ge 0.95$ (Figure 5); for greater $\overline{\mathrm{A}}$ values, the return periods are much smaller for conditional $\text{LD}$ compared to smaller $\overline{\mathrm{A}}$ values.

4. Conclusions

Droughts have many different spatial–temporal characteristics; hence, a better way to describe them is to use joint probability-based functions, rather than building drought-related analyses on univariate datasets as joint distribution of these drought characteristics often do not follow a particular known distribution. The current study offers a methodology to analyze the spatial–temporal multivariate aspect of droughts via copula functions in a probabilistic approach.

The results show that, among the drought characteristics, $\overline{\mathrm{S}}$ and $\overline{\mathrm{A}}$ are found to have significant linear dependence, implying their predictions using each other may yield reasonable forecasts, while the same may not be true for other relationships including $\text{LD}$. On the other hand, bivariate $\mathrm{T}$ is more sensitive to variability in $\text{LD}$ than $\overline{\mathrm{S}}$ or $\overline{\mathrm{A}}$, implying that predictions of bivariate $\mathrm{T}$ could be more accurate with the presence of $\text{LD}$ than $\overline{\mathrm{S}}$ or $\overline{\mathrm{A}}$. Nevertheless, observations of $\overline{\mathrm{S}}$ or $\overline{\mathrm{A}}$ narrow the variability window of trivariate $\mathrm{T}$ predictions, particularly for higher values of $\overline{\mathrm{S}}$ or $\overline{\mathrm{A}}$; implying that $\overline{\mathrm{S}}$ or $\overline{\mathrm{A}}$ observations have utility in predictions of $\mathrm{T}$ even though such predictions are more sensitive to variability in $\text{LD}$.

The flexibility of $\mathrm{T}$ calculation for different multivariate drought scenarios clearly illustrates the power of the introduced methodology. Additionally, certain drought characteristics may precede others, hence such observed drought characteristics may prove very valuable when making inferences about non-observed drought characteristics, while such relations can be easily studied using the multivariate conditional copula methodology presented in this study. In particular, trivariate drought analysis scenarios could prove very helpful in estimating a particular drought characteristic given availability of other drought characteristic info.

Drought characteristics often have dataset-specific distributions. Even though elliptical family-based copula functions yielded better performance in this study, other studies found that different copula families perform better [13,20]. This implies that the optimality of the copula functions could be dataset- and case-specific; hence inter-comparisons of copula functions could be necessary before estimation of return periods in other studies. However, this result should be confirmed via independent studies before it can be generalized for copula function-related drought analyses.

In this study, ground station-based datasets are utilized in all of the analyses. While they provide the longest datasets over the most locations, these station-based datasets are spatially limited, often contain gaps in their time series, and may have user-dependent data quality. On the other hand, remote sensing-based datasets prove very powerful alternatives to station-based datasets as these spatially and temporally consistent datasets are available everywhere, including remote locations where installation of stations is not economically feasible. Therefore, remotely sensed datasets could be recommended as alternative reanalysis products to be used as such model datasets may provide consistent and sufficiently long time series for drought analysis.