Copula-Based Severity–Duration–Frequency (SDF) Analysis of Streamﬂow Drought in the Source Area of the Yellow River, China

: In classical severity–duration–frequency (SDF) analysis, the dependence between drought characteristics is not effectively considered. The present study aims to propose the SDF relationships of streamﬂow drought in the source area of the Yellow River (SAYR) using a copula-based approach. Comparison of multiple time-varying threshold levels and the integration and elimination of drought events were considered. Selection of marginal probability distribution and copula-based joint probability distribution was properly conducted with multiple means. Copula-based joint and conditional probabilities were computed. The ﬁndings support carrying out integration and elimination processing on the preliminarily identiﬁed streamﬂow droughts through a run analysis with a time-varying threshold level of the 80% quantile of daily streamﬂow. The Gaussian copula was selected as the optimal model for constructing bivariate joint probability distribution, with generalized extreme value and log-normal as the suitable marginal probability distributions of streamﬂow drought duration and severity. The proposed copula-based SDF relationships of streamﬂow drought events can provide more critical information in addition to univariate frequency analysis, beneﬁtting from the joint and conditional probabilities. The multivariate probabilistic analyses can effectively consider the connection and interaction between drought characteristics, while conditional probability distribution allows analyzing the impact of one drought characteristic on another. The results also indicate a relatively high risk of streamﬂow drought with short duration and low severity in the region, requiring effective drought-mitigation strategies and measures.


Introduction
Drought is a kind of natural phenomenon with complicated causes, mainly manifested as abnormal water supply conditions. From different research perspectives, droughts are usually divided into meteorological, agricultural, hydrological, socio-economic, ecological, and groundwater categories [1,2]. Among them, hydrological drought refers to the reduction of surface water due to lack of precipitation, such as river flow and lake/reservoir storage, which may cause serious consequences [3]. Streamflow drought is the focus of hydrological drought research, and it is usually possible to identify streamflow drought events based on the theory of runs [4]. In run analysis, the truncation threshold can be either fixed or time-varying. However, considering the significant difference in river flow between wet and dry seasons, it is unreasonable to use a fixed threshold to identify streamflow drought during a long period [5]. At the same time, in identifying drought events, there are often cases where the time interval between adjacent drought events is very short as well

Study Area and Data
The source area of the Yellow River (SAYR) is located in the northeast of the Qinghai-Tibet Plateau. As the gateway of the SAYR, the Tangnaihai gauge (100 • 09 E, 35 • 30 N) is situated on the main stream of the Yellow River at Tangnaihai Township, Xinghai County, Qinghai Province, China, providing important basic data about variations of hydrological processes and water resources as well as changes in the water-sediment relationship in the upper reaches of the Yellow River. The geographical features of the SAYR and relevant catchment area drained at the Tangnaihai gauge are shown in Figure 1. In this study, the observational datasets of daily averaged streamflow at the Tangnaihai gauge from 1956 to 2018 were used to analyze the historical streamflow drought characteristics (e.g., duration, severity, and severity peak) and their joint probabilistic relationships from a multivariate perspective in the SAYR.

Study Area and Data
The source area of the Yellow River (SAYR) is located in the northeast of the Qinghai-Tibet Plateau. As the gateway of the SAYR, the Tangnaihai gauge (100°09′ E, 35°30′ N) is situated on the main stream of the Yellow River at Tangnaihai Township, Xinghai County, Qinghai Province, China, providing important basic data about variations of hydrological processes and water resources as well as changes in the water-sediment relationship in the upper reaches of the Yellow River. The geographical features of the SAYR and relevant catchment area drained at the Tangnaihai gauge are shown in Figure 1. In this study, the observational datasets of daily averaged streamflow at the Tangnaihai gauge from 1956 to 2018 were used to analyze the historical streamflow drought characteristics (e.g., duration, severity, and severity peak) and their joint probabilistic relationships from a multivariate perspective in the SAYR.

Streamflow Drought Identification
To start with, a series of abnormally low streamflow episodes can be identified based on the run theory and using a given flow rate as a threshold level. As a useful tool, the run analysis is commonly used to examine the internal coherence of a time-dependent sequence based on the arrangement of sample markers (i.e., adjacently identical or different). Through run analyses, a negative run process is defined as a streamflow drought when daily streamflow is consistently lower than the threshold level during a continuous period. As illustrated in Figure 2, a streamflow drought is generally characterized by duration ( , d), severity ( , 10 8 m 3 ), severity peak ( , m 3 /s), and inter-arrival time (d). Specifically, duration refers to the length of time from the beginning to the end of a streamflow drought, severity records the accumulative water shortage (in the form of volume) relevant to drought duration, severity peak reflects the maximum negative anomaly of streamflow during an identified drought episode, and inter-arrival time indicates the time span between two neighboring streamflow droughts. According to previous studies (e.g., [5]), the 70% to 95% quantiles of daily streamflow (i.e., Q70 to Q95) are usually used as candidates of threshold levels for run analyses in streamflow drought identification. In this way, the daily threshold becomes time-varying and is composed of 365 or 366 threshold values per year, which could be derived from daily streamflow datasets of a relatively long record through frequency calculation. Figure 3 shows the illustration of streamflow drought identification using a time-varying threshold level of the 80% quantile of daily

Streamflow Drought Identification
To start with, a series of abnormally low streamflow episodes can be identified based on the run theory and using a given flow rate as a threshold level. As a useful tool, the run analysis is commonly used to examine the internal coherence of a time-dependent sequence based on the arrangement of sample markers (i.e., adjacently identical or different). Through run analyses, a negative run process is defined as a streamflow drought when daily streamflow is consistently lower than the threshold level during a continuous period. As illustrated in Figure 2, a streamflow drought is generally characterized by duration (D, d), severity (S, 10 8 m 3 ), severity peak (P, m 3 /s), and inter-arrival time (d). Specifically, duration refers to the length of time from the beginning to the end of a streamflow drought, severity records the accumulative water shortage (in the form of volume) relevant to drought duration, severity peak reflects the maximum negative anomaly of streamflow during an identified drought episode, and inter-arrival time indicates the time span between two neighboring streamflow droughts. According to previous studies (e.g., [5]), the 70% to 95% quantiles of daily streamflow (i.e., Q70 to Q95) are usually used as candidates of threshold levels for run analyses in streamflow drought identification. In this way, the daily threshold becomes time-varying and is composed of 365 or 366 threshold values per year, which could be derived from daily streamflow datasets of a relatively long record through frequency calculation. Figure 3 shows the illustration of streamflow drought identification using a time-varying threshold level of the 80% quantile of daily streamflow (Q80). It is seen that the daily streamflow fluctuates below and above the threshold level from time to time during a calendar year, while a series of streamflow droughts and relevant characteristics (e.g., duration, severity, and severity peak) can be extracted regarding a daily streamflow lower than Q80. Meanwhile, more and relatively longer streamflow droughts are generally identified during low-flow season.
Water 2023, 15, x FOR PEER REVIEW 4 of 20 streamflow (Q80). It is seen that the daily streamflow fluctuates below and above the threshold level from time to time during a calendar year, while a series of streamflow droughts and relevant characteristics (e.g., duration, severity, and severity peak) can be extracted regarding a daily streamflow lower than Q80. Meanwhile, more and relatively longer streamflow droughts are generally identified during low-flow season.

Integration and Elimination of Drought Events
For daily streamflow series, there are often cases where the streamflow temporarily exceeds the threshold level (theoretically non-drought) between two identified adjacent streamflow droughts with relatively long duration and high severity. In this situation, it is usually considered that there might be a connection between these two artificially fragmented droughts, and they should be integrated into one drought event of a longer duration and higher magnitude (see Figure 4). Ref. [33] introduced a method for integrating drought events according to the interval and over-threshold flow rate between two adjacent droughts. Assume that there are two adjacent streamflow droughts, denoted by and , whose duration, severity, and severity peak are , , and , , , respectively, and they should be integrated if the following conditions are met: (1) where is the time interval between two adjacent droughts; is the critical value of , assumed to be 5 days in this study; is the sum of over-threshold flow; is the severity of streamflow drought ; is the critical value of (ratio of and ), generally taking a value of 0.1.

Integration and Elimination of Drought Events
For daily streamflow series, there are often cases where the streamflow temporarily exceeds the threshold level (theoretically non-drought) between two identified adjacent streamflow droughts with relatively long duration and high severity. In this situation, it is usually considered that there might be a connection between these two artificially fragmented droughts, and they should be integrated into one drought event of a longer duration and higher magnitude (see Figure 4). Ref. [33] introduced a method for integrating drought events according to the interval and over-threshold flow rate between two adjacent droughts. Assume that there are two adjacent streamflow droughts, denoted by R i and R i+1 , whose duration, severity, and severity peak are (D i , S i , P i ) and (D i+1 , S i+1 , P i+1 ), respectively, and they should be integrated if the following conditions are met: where T i is the time interval between two adjacent droughts; T c is the critical value of T i , assumed to be 5 days in this study; O i is the sum of over-threshold flow; S i is the severity of streamflow drought R i ; ρ c is the critical value of ρ i (ratio of O i and S i ), generally taking a value of 0.1. Supposing an integrated streamflow drought, , the drought integration is specifically conducted as follows: , where , , and are the duration, severity, and severity peak of the integrated streamflow drought , respectively. For integrated and following streamflow droughts, the criteria described in Equation (1) are examined as well. To be specific, if the interval and over-threshold flow rate have not reached given critical values, the integration of these two adjacent droughts are further carried out using the above procedure given by Equation (2). In other words, the integration of streamflow droughts would iterate until the requirements of the interval and over-threshold flow rate are no longer met. For example, the three streamflow droughts in the left part of Figure 4 (with severities of , , and , and sums of the over-threshold flow of and ) are integrated into only one drought event with a two-round integration, of which the duration, severity, and severity peak can be derived according to Equation (2).
After the integration of drought events, there might still be many drought events of short duration or low severity. These minor drought events are of little statistical significance for drought analysis and are suggested to be eliminated, such as the fourth drought event with the severity of in Figure 4. Assuming the mathematical expectations of duration and severity of all drought events are and , respectively, a streamflow drought event should be eliminated if its duration is shorter than , or its severity is smaller than , where 0.3 are predetermined coefficients [33].

Marginal Probability Distribution
The characteristics (duration, severity, and severity peak) of streamflow drought events after integration and elimination treatments, as random variables, are separately fitted to appropriate probability distributions that can produce good agreements between their empirical and theoretical frequencies. To this end, the L-moment ratio diagrams are used to select the most suitable distributions out of quite a few two-and three-parameter candidate ones. As a useful visual inspection tool, the L-moment ratio diagrams can provide comparative information on shape features (i.e., L-moments that specify a unique Supposing an integrated streamflow drought, R U , the drought integration is specifically conducted as follows: where D U , S U , and P U are the duration, severity, and severity peak of the integrated streamflow drought R U , respectively. For integrated and following streamflow droughts, the criteria described in Equation (1) are examined as well. To be specific, if the interval and over-threshold flow rate have not reached given critical values, the integration of these two adjacent droughts are further carried out using the above procedure given by Equation (2). In other words, the integration of streamflow droughts would iterate until the requirements of the interval and overthreshold flow rate are no longer met. For example, the three streamflow droughts in the left part of Figure 4 (with severities of S 1 , S 2 , and S 3 , and sums of the over-threshold flow of O 1 and O 2 ) are integrated into only one drought event with a two-round integration, of which the duration, severity, and severity peak can be derived according to Equation (2).
After the integration of drought events, there might still be many drought events of short duration or low severity. These minor drought events are of little statistical significance for drought analysis and are suggested to be eliminated, such as the fourth drought event with the severity of S 4 in Figure 4. Assuming the mathematical expectations of duration and severity of all drought events are E{ D} and E{ S} , respectively, a streamflow drought event should be eliminated if its duration is shorter than r d E{ D} , or its severity is smaller than r s E{ S} , where r d = r s = 0.3 are predetermined coefficients [33].

Marginal Probability Distribution
The characteristics (duration, severity, and severity peak) of streamflow drought events after integration and elimination treatments, as random variables, are separately fitted to appropriate probability distributions that can produce good agreements between their empirical and theoretical frequencies. To this end, the L-moment ratio diagrams are used to select the most suitable distributions out of quite a few two-and three-parameter candidate ones. As a useful visual inspection tool, the L-moment ratio diagrams can provide comparative information on shape features (i.e., L-moments that specify a unique distribution, especially the L-skewness τ 3 and L-kurtosis τ 4 ) of empirical and theoretical distributions [34,35]. The marginal probability distributions selected for streamflow drought duration, severity, and severity peak are given, respectively, as follows.
(1) The generalized extreme value (GEV) distribution with a probability density function: (2) The log-normal distribution with a probability density function: (3) The generalized Pareto distribution with a probability density function: where x stands for the drought variable (duration, severity, or severity peak); µ, σ, and κ are location, scale, and shape parameters of relevant probability distributions, respectively, which can be estimated through the maximum likelihood approach. The goodness-of-fit of streamflow drought characteristics (duration, severity, and severity peak) fitted with selected theoretical distributions are also examined by root mean square error (RMSE) and Kolmogorov-Smirnov (K-S) test. The RMSE is an effective bias statistic based on mean square error (MSE): where n is the sample size, and m is the number of parameters in the considered probability distribution or statistical model; p e and p are the probabilities calculated through the empirical plotting-position formula and derived from the selected theoretical model, respectively. In general, a small RMSE value indicates little fitting bias and good fitness between empirical points and theoretical model. As a nonparametric statistical method, the sample statistic K n of K-S test can be calculated as given: where n is the sample size, and p (•) stands for the theoretical probabilities based on an ascendingly sorted sample: Given a statistically significant level (e.g., α = 0.05), the critical statistic K 0 of the K-S test is obtained through Monte Carlo stochastic simulation, with a bootstrap repetition of 5000 times. If the sample statistic K n is less than the critical statistic K 0 , the null hypothesis of the K-S test is accepted. In other words, the goodness-of-fit of the observed sample fitted with the checked theoretical model is supported by the K-S test.

Joint Probability Distribution
Copulas are sophisticated tools connecting marginal probability distributions of dependent variables to produce multivariate joint probability distribution. The dependence structure of correlated variables is implicitly considered by appropriate copula functions. Thus, correlationship analysis is needed to preclude independent cases. Both linear (Pear- son's γ n ) and nonlinear rank-based (Spearman's ρ n and Kendall's τ n ) correlation coefficients are used to measure the paired dependence of streamflow drought characteristics (i.e., duration-severity, duration-peak, and severity-peak). A larger correlation coefficient suggests a relatively higher dependent relationship between analyzed variables. At the same time, the paired dependence of streamflow drought characteristics (e.g., duration and severity) is also inspected using graphic analysis tools of Chi-plot and K-plot. The Chi-plot and K-plot are both based on the rank of sample points and can diagnose potential relationships through comparing the distribution of sample points against the theoretical dependent area [36,37].
Two meta-elliptical (Gaussian and Student's t) and three Archimedean (Frank, Gumbel, and Clayton) copulas are used to construct a bivariate joint probability distribution of significantly dependent streamflow drought characteristics (e.g., duration and severity). Among them, the expression of Gaussian copula is given below, while readers are referred to other literature (e.g., [38,39]) for more details on other copula functions. Specifically, the cumulative distribution function of a two-dimensional Gaussian copula is as follows: where C(•) stands for the cumulative probability of bivariate joint distribution; u 1 = F D (d) and u 2 = F S (s) are the marginal probability distributions of streamflow drought duration (D) and severity (S), respectively; is the symmetrical covariance matrix, and ρ 12 = ρ 21 = ρ is the parameter of two-dimensional Gaussian copula to be estimated through the maximum likelihood approach; w = [w 1 , w 2 ] T represents the vector of integral variables. In addition to the RMSE, the Akaike information criterion (AIC) and Bayesian information criterion (BIC) are also adopted to assess the goodness-of-fit of constructed bivariate joint distribution using different copula functions. The AIC and BIC, taking into account both model-fitting bias and the number of parameters, are, respectively, formulated as follows: AIC = n ln(MSE) + 2m (10) where n is the sample size; m is the number of model parameters; 2m and m ln(n) are the penalty terms; MSE is as defined in Equation (6). The minimum value of AIC and/or BIC values indicates the best fitness between sample points and the theoretical model. Moreover, the K-S test with a modified bootstrap procedure based on Rosenblatt's transformation is also conducted to validate the goodness-of-fit of the selected copula function. Analogous to the K-S test of marginal probability distribution, the theoretical and empirical probabilities in Equation (8) are, respectively, derived from bivariate copula and plotting-position formula for joint probability distribution, whose difference provides a critical basis for calculating the statistics of K-S test. For a detailed methodology and procedure, one may refer to [28].

Copula-Based Joint and Conditional Probabilities
Since the joint distribution of streamflow drought characteristics (e.g., duration and severity) is finally described with selected copula function, the bivariate joint non-exceedance probability can be derived as given: where F D,S (d, s) is the joint probability distribution of drought duration and severity, with F D (d) and F S (s) being the marginal probability distribution functions; C(u 1 , u 2 ) denotes the cumulative distribution function of the two-dimensional copula, while u 1 = F D (d) and u 2 = F S (s). Accordingly, the bivariate joint exceedance probability is then obtained: Based on the bivariate joint probability distribution, the conditional non-exceedance probability of the streamflow drought duration (severity) given severity (duration) exceeding a certain magnitude can also be estimated:

Drought Identification
In order to ascertain the impact of different threshold levels on identified streamflow droughts and determine the most appropriate threshold for the Tangnaihai gauge, six comparable thresholds were considered, namely Q70, Q75, Q80, Q85, Q90, and Q95 of daily streamflow, to extract the drought events together with corresponding duration (D), severity (S), and severity peak (P), whose statistics are compared in Table 1. It was found that with the gradual reduction of daily thresholds from Q70 to Q95, the number of droughts (including short-duration droughts) and the averages of drought duration and severity generally show a downward trend. However, if the threshold is too small (e.g., Q95), although the proportion of short-duration droughts is relatively higher, the averages of drought duration and severity are larger than those of the daily threshold of Q90. It implies that a very small threshold level would lead to a few drought events with considerably long duration and huge severity, thereby failing to reflect the essential statistics of streamflow drought characteristics. Meanwhile, although the number of shortduration droughts declines with decreasing threshold levels, its proportion overall increases. In other words, the use of too-small threshold levels would result in a high fraction of short-duration droughts having little significance but causing interference with statistical analysis of streamflow droughts. Therefore, too-small threshold levels should be avoided for streamflow drought identification. It was also seen that the daily threshold level of Q80 yields the smallest proportion of short-duration droughts. Thus, the following parts mainly focus on analyses of streamflow droughts at the Tangnaihai gauge using Q80 of daily streamflow as a selected time-varying threshold level. Drought integration and elimination were further conducted on preliminarily identified streamflow drought events with a selected daily threshold level of Q80. The statistics of streamflow drought events with different treatments (untreated, only integration, and integration and elimination) are displayed in Table 2. It can be seen that with no treatment, up to 429 drought events were identified at the Tangnaihai gauge from 1956 to 2018, i.e., the occurrence of a drought event every 1.76 months on average. After integration and elimination, a total of 93 drought events were concentrated instead, i.e., a drought event occurring about every 8.13 months. Accordingly, the average of drought duration increases from 10.3 to 41.2 days. It is considered that the statistics of duration and severity with drought integration and elimination can better reflect the properties of streamflow drought occurrence and persistence. Therefore, it is of great significance to integrate and eliminate identified streamflow drought events before fitting the drought characteristics (e.g., duration, severity, and severity peak) to suitable probability distributions and analyzing their joint probabilistic relationships.

Selection of Marginal Probability Distributions
It is convenient and visible to select appropriate probability distributions for streamflow drought duration, severity, and severity peak using the L-moment ratio diagram (see Figure 5). To be specific, the ratio of given L-moments (i.e., L-skewness τ 3 and L-kurtosis τ 4 ) of normal, exponential, and Gumbel distributions as points with fixed positions and generalized Pareto, log-logistic, GEV, Pearson type III, and log-normal distributions as a cluster of curves were drawn as reference. Then, the sample L-moments (L-skewness τ 3 and L-kurtosis τ 4 ) of drought duration, severity, and severity peak were calculated and displayed against the theoretical L-moments of the above candidate probability distributions. As shown in Figure 5, the ratio of L-skewness τ 3 and L-kurtosis τ 4 of drought duration (D) is close to the theoretical curve of GEV or log-logistic distribution, while log-normal and generalized Pareto distributions stand out in generating the sample L-moments of drought severity (S) and severity peak (P), respectively. By contrast, the theoretical L-moments of other probability distributions are far away from the sample L-moments of drought duration (D) and severity (S), such as normal, exponential, and Gumbel distributions.
The cumulative probabilities of streamflow drought duration and severity based on GEV and log-normal distributions are demonstrated in Figure 6. As can be seen, the theoretical probabilistic plots of GEV and log-normal distributions generally produce an acceptable fitness to the empirical frequency of drought duration and severity, respectively. Moreover, the theoretical and empirical probabilities are also compared using the P-P plots, which generally show an agreement with most data points near the 45-degree line (especially the upper tail, with a larger duration or severity that is of greater significance for drought analysis).
For streamflow drought characteristics fitted by different probability distributions (i.e., GEV for duration, log-normal for severity, and generalized Pareto for severity peak), the distribution parameters are estimated using the maximum likelihood approach and given in Table 3. By doing so, the probabilistic properties of drought duration, severity, and severity peak are largely characterized by selected probability distributions with estimated location, scale, and shape parameters. Meanwhile, the RMSE values, showing relative biases between theoretical and empirical probabilities, are very small and suggest a good fitness of drought characteristics to corresponding probability distributions. Furthermore, the selection of marginal probability distributions for drought duration, severity, and severity peak is also supported by the Kolmogorov-Smirnov (K-S) goodness-of-fit tests since all sample statistics are smaller than relevant critical values at the significant level of 0.05 (see the last two columns of Table 3). The cumulative probabilities of streamflow drought duration and severity based on GEV and log-normal distributions are demonstrated in Figure 6. As can be seen, the theoretical probabilistic plots of GEV and log-normal distributions generally produce an acceptable fitness to the empirical frequency of drought duration and severity, respectively. Moreover, the theoretical and empirical probabilities are also compared using the P-P plots, which generally show an agreement with most data points near the 45-degree line (especially the upper tail, with a larger duration or severity that is of greater significance for drought analysis). For streamflow drought characteristics fitted by different probability distributions (i.e., GEV for duration, log-normal for severity, and generalized Pareto for severity peak), the distribution parameters are estimated using the maximum likelihood approach and given in Table 3. By doing so, the probabilistic properties of drought duration, severity, and severity peak are largely characterized by selected probability distributions with estimated location, scale, and shape parameters. Meanwhile, the RMSE values, showing relative biases between theoretical and empirical probabilities, are very small and suggest a good fitness of drought characteristics to corresponding probability distributions. Furthermore, the selection of marginal probability distributions for drought duration, severity, and severity peak is also supported by the Kolmogorov-Smirnov (K-S) goodness-offit tests since all sample statistics are smaller than relevant critical values at the significant level of 0.05 (see the last two columns of Table 3).

Construction of Bivariate Joint Probability Distribution
Cross-correlation analysis between variables is required before using the copula function to construct the joint probability distribution of drought characteristics. The results of multiple correlation coefficients, including Pearson's γ n , Spearman's ρ n , and Kendall's τ n , were computed and are presented in Table 4. It is indicated that drought duration (D) and severity (S) are highly correlated, with Pearson's linear correlation coefficient γ n being up to 0.85. However, drought severity (S) and severity peak (P) seem to have a relatively higher correlation from a nonlinear perspective, with larger correlation coefficients of Spearman's ρ n . On the other hand, drought duration (D) and severity peak (P) show a relatively low correlation in terms of both linear and nonlinear correlation coefficients. Table 4. Paired correlation coefficients of streamflow drought characteristics. In addition, the significant correlationship between drought duration (D) and severity (S) is also visible in the rank-based Chi-plot and K-plot (see Figure 7). In the Chi-plot, the area between two dotted horizontal lines of ±0.18 is defined as the confidence band, with a significant level of 0.05, which means an independent relationship. Since most of the data points fall outside the confidence band in Chi-plot, it can be generally concluded that apparent dependence exists between drought duration (D) and severity (S). Meanwhile, the dotted vertical zero line and the distribution of data points also imply a positive asymmetric correlation structure between the two drought characteristics. On the contrary, the range enclosed by the dotted curve and diagonal line in the K-plot indicates a dependent relationship. Further, the closer the data points to the curve, the higher the suggested dependence. It is seen that almost all the data points fall inside the dependent area and closer to the indicative curve for the upper tail. Overall, the correlation and dependence between drought duration (D) and severity (S) are validated, and thereby, it is possible to construct their joint probability distribution using appropriate copula functions. coefficients of Spearman's . On the other hand, drought duration ( ) and severity peak ( ) show a relatively low correlation in terms of both linear and nonlinear correlation coefficients. In addition, the significant correlationship between drought duration ( ) and severity ( ) is also visible in the rank-based Chi-plot and K-plot (see Figure 7). In the Chi-plot, the area between two dotted horizontal lines of 0.18 is defined as the confidence band, with a significant level of 0.05, which means an independent relationship. Since most of the data points fall outside the confidence band in Chi-plot, it can be generally concluded that apparent dependence exists between drought duration ( ) and severity ( ). Meanwhile, the dotted vertical zero line and the distribution of data points also imply a positive asymmetric correlation structure between the two drought characteristics. On the contrary, the range enclosed by the dotted curve and diagonal line in the K-plot indicates a dependent relationship. Further, the closer the data points to the curve, the higher the suggested dependence. It is seen that almost all the data points fall inside the dependent area and closer to the indicative curve for the upper tail. Overall, the correlation and dependence between drought duration ( ) and severity ( ) are validated, and thereby, it is possible to construct their joint probability distribution using appropriate copula functions. For the construction of bivariate joint distribution of drought duration and severity, the maximum likelihood estimates of parameters of five two-dimensional copulas (i.e., Gaussian, Student's t, Frank, Gumbel, and Clayton) are provided in Table 5. Using these estimated parameters, the theoretical bivariate joint probabilities can be derived from the cumulative distribution functions of different copulas (e.g., Gaussian). Comparing to bivariate empirical probabilities based on the plotting-position formula, the fitting biases in terms of RMSE were calculated. Moreover, AIC and BIC were further computed to compare the goodness-of-fit of different copulas as a joint probability distribution of drought duration and severity. As shown in Table 5, Gaussian copula generally has the best performance, resulting in the smallest RMSE, AIC, and BIC values. The P-P plot of the joint probability distribution of streamflow drought duration and severity using a Gaussian copula is illustrated in Figure 8 For the construction of bivariate joint distribution of drought duration and severity, the maximum likelihood estimates of parameters of five two-dimensional copulas (i.e., Gaussian, Student's t, Frank, Gumbel, and Clayton) are provided in Table 5. Using these estimated parameters, the theoretical bivariate joint probabilities can be derived from the cumulative distribution functions of different copulas (e.g., Gaussian). Comparing to bivariate empirical probabilities based on the plotting-position formula, the fitting biases in terms of RMSE were calculated. Moreover, AIC and BIC were further computed to compare the goodness-of-fit of different copulas as a joint probability distribution of drought duration and severity. As shown in Table 5, Gaussian copula generally has the best performance, resulting in the smallest RMSE, AIC, and BIC values. The P-P plot of the joint probability distribution of streamflow drought duration and severity using a Gaussian copula is illustrated in Figure 8, where horizontal and vertical axes show the empirical and theoretical probabilities derived from the plotting-position formula and bivariate copula (using estimated parameters presented in Table 5), respectively. It is seen that the data points are all along the 45-degree line, suggesting a good agreement between the bivariate empirical and theoretical probabilities. Moreover, the competence of Gaussian copula is still supported by the K-S goodness-of-fit test. Specifically, the sample statistic (0.0684) is less than the relevant critical value (0.1497), with a significant level of 0.05. Table 5. Comparison of two-dimensional copulas for constructing joint probability distribution of streamflow drought duration and severity.

Copulas
Parameters Parameter Estimates and theoretical probabilities derived from the plotting-position formula and bivariate copula (using estimated parameters presented in Table 5), respectively. It is seen that the data points are all along the 45-degree line, suggesting a good agreement between the bivariate empirical and theoretical probabilities. Moreover, the competence of Gaussian copula is still supported by the K-S goodness-of-fit test. Specifically, the sample statistic (0.0684) is less than the relevant critical value (0.1497), with a significant level of 0.05. Note: * The degree of freedom for Student's t copula is 6. Figure 8. The P-P plot of joint probability distribution of streamflow drought duration and severity using Gaussian copula.

Gaussian Copula-Based SDF Relationships of Streamflow Drought
The severity-duration-frequency (SDF) relationships of streamflow drought in the SAYR were finally developed based on the two-dimensional Gaussian copula. For this purpose, the bivariate joint non-exceedance and exceedance probability as well as relevant conditional non-exceedance probability were computed on the basis of Equations (12)- (15). The probabilistic relationship among streamflow drought duration, severity, and frequency (in terms of bivariate joint probability) is illustrated by a three-dimensional surface and corresponding contours (see Figure 9). Accordingly, it is possible to obtain the joint probability with different combinations of drought characteristics. For instance, a combination of 100-day duration and 8 × 10 8 m 3 severity corresponds to a non-exceedance probability of about 0.9, while many other combinations of different duration and severity values would also produce the same non-exceedance probability for streamflow drought analysis. It is very different from drought frequency analysis, which considers only one variable (e.g., duration or severity) and fails to consider the connection and interaction between multiple drought characteristics well. Specifically, the univariate-based design Theoretical probability Figure 8. The P-P plot of joint probability distribution of streamflow drought duration and severity using Gaussian copula.

Gaussian Copula-Based SDF Relationships of Streamflow Drought
The severity-duration-frequency (SDF) relationships of streamflow drought in the SAYR were finally developed based on the two-dimensional Gaussian copula. For this purpose, the bivariate joint non-exceedance and exceedance probability as well as relevant conditional non-exceedance probability were computed on the basis of Equations (12)- (15). The probabilistic relationship among streamflow drought duration, severity, and frequency (in terms of bivariate joint probability) is illustrated by a three-dimensional surface and corresponding contours (see Figure 9). Accordingly, it is possible to obtain the joint probability with different combinations of drought characteristics. For instance, a combination of 100-day duration and 8 × 10 8 m 3 severity corresponds to a non-exceedance probability of about 0.9, while many other combinations of different duration and severity values would also produce the same non-exceedance probability for streamflow drought analysis. It is very different from drought frequency analysis, which considers only one variable (e.g., duration or severity) and fails to consider the connection and interaction between multiple drought characteristics well. Specifically, the univariate-based design value is just equivalent to the intersection of concerned probability contour and horizontal/vertical axis. For example, the design severity corresponding to a non-exceedance probability of 0.9 turns out to be about 6 × 10 8 m 3 using marginal probability distribution of duration in traditional SDF analysis of droughts. Likewise, the design duration of the same nonexceedance probability is about 80 days according to marginal probability distribution of severity. In other words, with an identical probability level, the traditional SDF curves would significantly underestimate the design values of drought characteristics (e.g., duration and severity) relative to the copula-based SDF analysis that is based on joint rather than marginal probability distribution. Similarly, the concurrent drought risk of both duration and severity exceeding specific values can also be estimated by the Gaussian copula-based bivariate exceedance probability and its contours, as shown in Figure 9. For example, the exceedance probability of drought duration and severity simultaneously exceeding 60 days and 3.5 × 10 8 m 3 , respectively, is about 0.1. In general, the bivariate exceedance probability when drought duration and severity reach certain levels (e.g., 100-day duration and 8.5 × 10 8 m 3 severity, respectively) becomes relatively small, such as less than 0.05. The bivariate joint exceedance probability focuses on the concurrent risk of both drought duration and severity exceeding critical alert levels, which might trigger a failure in water supply and/or destruction of infrastructure. On the other hand, the joint non-exceedance probability indicates the possibility of both drought duration and severity within tolerable limits, which provides an implication about the conditions of safety. A combined interpretation of them can help to achieve a more comprehensive evaluation of exposed risks to drought. value is just equivalent to the intersection of concerned probability contour and horizontal/vertical axis. For example, the design severity corresponding to a non-exceedance probability of 0.9 turns out to be about 6 × 10 8 m 3 using marginal probability distribution of duration in traditional SDF analysis of droughts. Likewise, the design duration of the same non-exceedance probability is about 80 days according to marginal probability distribution of severity. In other words, with an identical probability level, the traditional SDF curves would significantly underestimate the design values of drought characteristics (e.g., duration and severity) relative to the copula-based SDF analysis that is based on joint rather than marginal probability distribution. Similarly, the concurrent drought risk of both duration and severity exceeding specific values can also be estimated by the Gaussian copula-based bivariate exceedance probability and its contours, as shown in Figure 9. For example, the exceedance probability of drought duration and severity simultaneously exceeding 60 days and 3.5 × 10 8 m 3 , respectively, is about 0.1. In general, the bivariate exceedance probability when drought duration and severity reach certain levels (e.g., 100day duration and 8.5 × 10 8 m 3 severity, respectively) becomes relatively small, such as less than 0.05. The bivariate joint exceedance probability focuses on the concurrent risk of both drought duration and severity exceeding critical alert levels, which might trigger a failure in water supply and/or destruction of infrastructure. On the other hand, the joint nonexceedance probability indicates the possibility of both drought duration and severity within tolerable limits, which provides an implication about the conditions of safety. A combined interpretation of them can help to achieve a more comprehensive evaluation of exposed risks to drought. In addition, Figure 10 provides the conditional non-exceedance probability of the drought duration (severity) given severity (duration) exceeding a series of relevant levels. As shown, the conditional non-exceedance probability of drought duration and severity decreases with the increase in the corresponding conditional factor. It is also valuable and interesting to compare the marginal, joint, and conditional non-exceedance probabilities  Figure 9. The Gaussian copula-based joint non-exceedance (left) and exceedance (right) probabilities and contours of streamflow drought duration and severity.
In addition, Figure 10 provides the conditional non-exceedance probability of the drought duration (severity) given severity (duration) exceeding a series of relevant levels. As shown, the conditional non-exceedance probability of drought duration and severity decreases with the increase in the corresponding conditional factor. It is also valuable and interesting to compare the marginal, joint, and conditional non-exceedance probabilities of bivariate streamflow drought characteristics. For example, Figure 11 details the case of drought duration and severity being 30 days and 2 × 10 8 m 3 , respectively. It is seen that a non-exceedance probability of about 0.6 is recognized through univariate frequency analysis with the given drought duration and severity. However, the conditional non-exceedance probability of either of them with the other one exceeding the given value reduces dramatically to about 0.3, whereas the bivariate joint non-exceedance probability of at least one of them not exceeding the given value(s) turns out to be around 0.45. The above-derived bivariate joint and conditional probabilities of drought characteristics can provide useful information for water resources planning and management, which would help to comprehensively evaluate the risk of malfunction for a specific streamflow drought event and further estimate the return period of a specific situation in which a water supply system fails to meet the requirement of water demand.
of bivariate streamflow drought characteristics. For example, Figure 11 details the case of drought duration and severity being 30 days and 2 × 10 8 m 3 , respectively. It is seen that a non-exceedance probability of about 0.6 is recognized through univariate frequency analysis with the given drought duration and severity. However, the conditional non-exceedance probability of either of them with the other one exceeding the given value reduces dramatically to about 0.3, whereas the bivariate joint non-exceedance probability of at least one of them not exceeding the given value(s) turns out to be around 0.45. The abovederived bivariate joint and conditional probabilities of drought characteristics can provide useful information for water resources planning and management, which would help to comprehensively evaluate the risk of malfunction for a specific streamflow drought event and further estimate the return period of a specific situation in which a water supply system fails to meet the requirement of water demand.

Discussion and Conclusions
For classical severity-duration-frequency (SDF) curves, the correlation between drought duration and severity is not effectively considered, which only provides, in essence, the severity-frequency relationship with given discrete values of duration or the duration-frequency relationship with given discrete values of severity. By contrast, the proposed copula-based SDF relationships take into account the dependence (e.g., degree and structure) between drought characteristics, providing superior severity-duration-frequency relationships that are more comprehensive and conform to reality. Using higherdimensional (e.g., three-dimensional) copula functions, more complicated multivariate relationships and curves of drought characteristics, such as four-dimensional duration-severity-peak-frequency relationships of streamflow drought events, could further be considered and constructed according to the framework and procedure adopted by the present study.
Using daily streamflow observations at the Tangnaihai gauge in the SAYR from 1956 to 2018, this study first obtained the time series of drought duration (D), severity (S), and severity peak (P) based on run analysis with a time-varying daily threshold and integration and elimination of drought events. Then, the bivariate joint probability distribution of drought duration and severity was constructed using two-dimensional copula func-

Discussion and Conclusions
For classical severity-duration-frequency (SDF) curves, the correlation between drought duration and severity is not effectively considered, which only provides, in essence, the severity-frequency relationship with given discrete values of duration or the duration-frequency relationship with given discrete values of severity. By contrast, the proposed copula-based SDF relationships take into account the dependence (e.g., degree and structure) between drought characteristics, providing superior severityduration-frequency relationships that are more comprehensive and conform to reality. Using higher-dimensional (e.g., three-dimensional) copula functions, more complicated multivariate relationships and curves of drought characteristics, such as fourdimensional duration-severity-peak-frequency relationships of streamflow drought events, could further be considered and constructed according to the framework and procedure adopted by the present study.
Using daily streamflow observations at the Tangnaihai gauge in the SAYR from 1956 to 2018, this study first obtained the time series of drought duration (D), severity (S), and severity peak (P) based on run analysis with a time-varying daily threshold and integration and elimination of drought events. Then, the bivariate joint probability distribution of drought duration and severity was constructed using two-dimensional copula functions with selected marginal probability distributions. Finally, the severity-duration-frequency (SDF) relationships based on Gaussian copula were developed and analyzed, mainly in terms of bivariate joint and conditional non-exceedance/exceedance probabilities. The main findings and conclusions are summarized as follows.
(1) The proportion of short-duration droughts generally increases as the threshold decreases, which suggests avoiding too-small thresholds, and the time-varying daily threshold level of Q80 is recommended for streamflow drought identification in the SAYR. After integration and elimination, the streamflow drought events are more consistent with the drought occurrence and persistence feature, highlighting the necessity to carry out integration and elimination processing on preliminarily identified streamflow droughts through run analysis; (2) According to the L-moment ratio diagram, P-P plot, and RMSE, the generalized extreme value, log-normal, and generalized Pareto are, respectively, suitable as marginal probability distributions of streamflow drought duration, severity, and severity peak at the Tangnaihai gauge. The correlation coefficient and rank-based correlation diagram suggested a significant asymmetric, positive correlation between drought duration and severity. Then, supported by the RMSE, AIC, and BIC, the Gaussian copula was selected as the optimal model for constructing bivariate joint probability distribution of streamflow drought duration and severity. In addition, the marginal and joint probability distributions of drought characteristics passed the K-S goodness-of-fit tests at the significant level of 0.05; (3) Compared to traditional SDF analysis, the proposed copula-based SDF relationships of streamflow drought events can provide more critical information. Specifically, with given non-exceedance/exceedance probabilities, it can consider the different combinations of multiple drought characteristics, making up for the defect of ignoring their connection and interaction in univariate frequency analysis. Thus, the corresponding multivariate probabilistic analyses are more comprehensive and more consistent with the essential attributes of drought events. Moreover, the conditional probability distribution effectively reflects the trend of gradually decreasing non-exceedance probabilities of drought duration (severity) with increasing severity (duration), which has practical significance for analyzing the probabilistic impact of one drought characteristic on another. From a multivariate perspective, the probability of one or two drought characteristics exceeding specific values would increase with decreasing non-exceedance but increasing exceedance probabilities. That is, the expected interarrival time of the designed drought event would be shorter as well. The results also indicate that the overall risk of streamflow drought with short duration and low severity is relatively high in the SAYR, and more attention is needed regarding effective drought-mitigation strategies and measures. Data Availability Statement: All data generated or analyzed are included in this published article.