Bivariate Assessment of Drought Return Periods and Frequency in Brazilian Northeast Using Joint Distribution by Copula Method

The Northeast region of Brazil (NRB) is the most populous semiarid area in the world and is extremely susceptible to droughts. The severity and duration of these droughts depend on several factors, and they do not necessarily follow the same behavior. The aim of this work is to evaluate the frequency of droughts in the NRB and calculate the return period of each drought event using the copula technique, which integrates the duration and severity of the drought in the NRB in a joint bivariate distribution. Monthly precipitation data from 96 meteorological stations spatially distributed in the NRB, ranging from 1961 to 2017, are used. The copula technique is applied to the Standardized Precipitation Index (SPI) on the three-month time scale, testing three families of Archimedean copula functions (Gumbel–Hougaard, Clayton and Frank) to reveal which model is best suited for the data. Averagely, the most frequent droughts observed in the NRB are concentrated in the northern sector of the region, with an observed duration varying from three and a half to five and a half months. However, the eastern NRB experiences the most severe droughts, lasting for 14 to 24 months. The probability distributions that perform better in modeling the series of severity and duration of droughts are exponential, normal and lognormal. The observed severity and duration values show that, for average values, the return period across the region is approximately 24 months. Still in this regard, the southernmost tip of the NRB stands out for having a return period of over 35 months. Regarding maximum observed values of severity and duration, the NRB eastern strip has the longest return period (>60 months), mainly in the southeastern portion where a return period above 90 months was observed. The northern NRB shows the shortest return period (~45 months), indicating that it is the NRB sector with the highest frequency of intense droughts. These results provide useful information for drought risk management in the NRB.


Introduction
Among adverse natural phenomena, drought is what most affects society by acting on large territorial extensions and lasting for long periods [1]. Overall, drought is characterized as a sufficiently long period of water deficiency to significantly impact edaphic, meteorological, hydrological and social aspects [2][3][4][5][6][7][8]. It is a phenomenon that imposes difficulties in the characterization of properties such as beginning, spatiotemporal extension and end, which makes its monitoring challenging [9].

Standardized Precipitation Index
Several indexes were created to study drought quantitatively. The Palmer Drought Severity Index (PDSI) was the first index developed for this purpose. Furthermore, PDSI was a very important index for incorporating the effects of evapotranspiration [52,53]. Since PDSI parameters were fixed We used a monthly accumulated precipitation time series from 96 conventional surface weather stations of the Brazilian National Institute of Meteorology (INMET), comprising the period from 1961 to 2017. The location of surface weather stations can be seen merged in Figure 1 (left panel). INMET data are subjected to a rigorous quality control and homogenization, providing confidence in data processing [47][48][49][50][51].
The concern with the quality of the rainfall data of this research is due to the fact that such observed series of precipitation will be used for validation of products from simulations of weather and climate associated with the project "Development and Validation of Modeling of the Components of the Brazilian Earth System Model (BESM)". This project has, among other objectives, the one of understanding the formation of clouds and precipitation which constitute one of the great challenges for the modeling of the terrestrial system and, consequently, for the studies of climate simulations.

Standardized Precipitation Index
Several indexes were created to study drought quantitatively. The Palmer Drought Severity Index (PDSI) was the first index developed for this purpose. Furthermore, PDSI was a very important index for incorporating the effects of evapotranspiration [52,53]. Since PDSI parameters were fixed and originally calculated for the USA, the self-calibrated PDSI (Sc-PDSI) was later defined, which estimated its parameters based on inherent series characteristics [54]. However, PDSI and Sc-PDSI did not allow analysis at various time scales. To cope with the problem of temporal characterization of droughts, the SPI was developed [55]. The calculation of SPI is based on long-term precipitation series at different time scales, such as 1, 3, 6, or 12 months. Such a property of the index makes it suitable to classify droughts that range from the 1-month edaphic type all the way up to hydrological ones that may last for years. The main premise of SPI is to assume that precipitation is the key variable for drought classification. More recently, a new index was developed to combine the main characteristics of SPI and PDSI, the Standard Precipitation and Evapotranspiration Index [9]. However, due to the lack of directly measured Evapotranspiration (Etp) data and the scarce availability of the data necessary to at least estimate it, we opted for the use of the SPI in the present study. Since its inception, several works have used the SPI for spatial and temporal drought studies [56][57][58], and several environmental centers have made the SPI the official index for monitoring the phenomenon. Moreover, the SPI is the drought index recommended by the World Meteorological Organization [59].
Given its widespread use and extensive literature validating the SPI, we chose to use the three-month time scale index (SPI-3) to derive Duration (D) and Severity (S) series. The D and S series are obtained from the Run Theory drought frequency analysis method ( Figure 2). Run Theory is a method for grouping independent drought events from time series based on the proximity of individual SPI observations with values below a threshold (0 in this work). Run Theory's premise is that very close negative SPI observations (even separated by small positive SPI intervals) represent observations from the same drought event. Added to the D and S of the events, the method provides the periods of non-occurrence of drought (nD) and the length (L) of periods of drought and non-drought. The mean values of nD and L will be used in Equation 3. Detailed information on SPI calculation is found in Mckee et al, [55]. Detailed information about the Run Theory can be found in Sun et al, [60].

Bivariate Copula
Copula functions are used to generate multivariate distributions from univariate marginal distributions. The theoretical basis of the technique was defined by Sklar. [61]. Subsequently, Genest et al., [62] established the definition of Archimedean copula functions, which often are used in hydrological variable analysis [32,[63][64][65]. observations from the same drought event. Added to the D and S of the events, the method provides the periods of non-occurrence of drought (nD) and the length (L) of periods of drought and nondrought. The mean values of nD and L will be used in Equation 3. Detailed information on SPI calculation is found in Mckee et.al, [55]. Detailed information about the Run Theory can be found in Sun et.al, [60].

2.3.Bivariate Copula
Copula functions are used to generate multivariate distributions from univariate marginal distributions. The theoretical basis of the technique was defined by Sklar. [61]. Subsequently, Genest et.al, [62] established the definition of Archimedean copula functions, which often are used in hydrological variable analysis [32,63−65]. Among the several existing Archimedean copula families, these work-tested copula functions of three families are frequently used in drought analyses (Table 1) [32,66,67]. Given that the series of S and D have a correlation (linear or nonlinear) to each other expressed by the Kendall correlation ) (τ , it is possible to estimate the parameter θ of the copula functions by the relationships shown in the third column of Table 1. We apply the Akaike Information Criterion (AIC) (Equation 2) to select, for each season, the copula that is best suited for modelling the properties of the drought [68]: Copula functions are defined as multivariate distribution functions, P(x 1 , . . . , x n ), with uniformly distributed variables on the interval [0,1]: where C is the Cumulative Distribution Function (CDF) of the copula, and P(x i ) is the marginal distribution of x i , being uniform on the interval [0,1], which also is denoted by u i The P(. . .) in the left-hand side is the joint probability distribution. Among the several existing Archimedean copula families, these work-tested copula functions of three families are frequently used in drought analyses (Table 1) [32,66,67]. Given that the series of S and D have a correlation (linear or nonlinear) to each other expressed by the Kendall correlation (τ), it is possible to estimate the parameter θ of the copula functions by the relationships shown in the third column of Table 1. We apply the Akaike Information Criterion (AIC) (Equation 2) to select, for each season, the copula that is best suited for modelling the properties of the drought [68]: where k is the number of the estimated parameter in the model andL is the maximum value of the likelihood function for the model. The copula function that has the lowest AIC value is the most suitable copula for the surface weather station. Table 1.
Archimedean copulas: families of tested functions (left column), functions of each family (central column) and relationship equations between τ and θ (right column).

Family Bivariate Copula Relationship between τ and θ
Gumbel-

Marginal Distribuitions for Drought Severity and Duration
A good performance in modeling the drought characteristics by joint distributions depends, fundamentally, on the appropriate choice of the marginal distributions that best fit the S and D data. Here, we tested the suitability of six probability distributions: normal, logistic, Weibull, gamma, lognormal and exponential. The choice of marginal distributions that best fit the data is based on the calculation of the Maximum Absolute Log-likelihood Value (MLV). Concerning each series of S and D, the MLV is calculated for the six distributions mentioned. The distribution showing the highest MLV is considered the most suitable for the series under consideration at the time.

Return Period Estimation
Drought return period analysis based on a single variable often leads to inaccurate results regarding event evolution. To model the drought return period more accurately, we performed our analysis based on bivariate joint distributions of duration and severity. The return period based on bivariate joint distributions can be defined as: where T (D≥dorS≥s) is the return period when one or both thresholds are known. E(L) is the expectation of the drought duration, which, in turn, is the sum of average D and nD ( Figure 2). Figure 3 shows the spatial distribution of mean and maximum drought severity over the NRB. The average drought severity across the NRB ranges from three to seven. The northern region has the highest values ( Figure 3a). Conversely, on average, the lowest values are observed in the east-south portion of the NRB. However, when it comes to the maximum severity values observed in the analysis period, the east-south portion and the entire east coast region showed the most severe droughts ( Figure 3b). Some points on the northern coast also exhibited extreme values relative to most of the NRB.

Drought Frequency
The northern NRB also stands out for its longer average drought duration, as D and S are positively correlated (Figure 4a). Further to the north of the NRB, average droughts of four and one half to seven months are noticed. The east-south portion of NRB, in turn, has the lowest D values, similar to the pattern observed for S (Figure 3a), with droughts lasting for two to three months. It is important to note that, although the two variables are highly correlated, the pattern will not necessarily be the same. Seen in northern NRB, a more uniform distribution of mean Severity is observed (Figure 3a), while mean Duration in the same sector exhibits a westward gradient (Figure 4a). Analyzing the maximum D observed throughout the NRB (Figure 4b), a pattern similar to maximum S (Figure 3b) is observed. The eastern region, which experienced the most severe droughts observed in the period under analysis, had the maximum Duration ever recorded, ranging from 14 to 24 months. The maximum Duration in the semiarid region varied from eight to 12 months.
The observed patterns are directly related to atmospheric phenomena that modulate climate variability in the NRB. The ITCZ is the main modulator of the rainy season in northern NRB [15,69,70]. During the austral summer, the ITCZ moves toward the southern hemisphere, due to the temperature gradient between the hemispheres, and reaches the northern portion of the NRB. A more southward displacement of the ITCZ causes increased precipitation, while a less southward migration of the same system causes a deficit of precipitation in the region. The main precipitation-generating phenomenon in the eastern portion of the NRB are the Easterly Wave Disturbances (EWD). This phenomenon consists of cloud bands that propagate westward in the trail of the trade winds and cause precipitation in the austral winter [45,46,71,72]. The frequency of this phenomenon determines the quality of the rainy season in the eastern sector. southward displacement of the ITCZ causes increased precipitation, while a less southward migration of the same system causes a deficit of precipitation in the region. The main precipitationgenerating phenomenon in the eastern portion of the NRB are the Easterly Wave Disturbances (EWD). This phenomenon consists of cloud bands that propagate westward in the trail of the trade winds and cause precipitation in the austral winter [45,46,71,72]. The frequency of this phenomenon determines the quality of the rainy season in the eastern sector.

Marginal Distributions
The multivariate distributions obtained by copula are generated from the marginal distributions of the individual variables. The advantage of the copula method is the flexibility to consider the probabilistic nature of each variable individually. Therefore, the marginal distributions are the probability distributions that best fit the variables under analysis. We tested the following probability distributions: normal (Gaussian), lognormal, exponential, logistic, Weibull and gamma [73]. According to the criterion of lowest AIC, the distributions that best fit the data are illustrated in Figure  5. The distributions selected for the severity series were the normal, lognormal and exponential ones (Figure 5a). The spatial pattern presented is similar to that observed in the fields of mean S and maximum S. Regarding the southeastern portion of the NRB, which had showed the lowest severity values (Figure 3a and Figure 4a), the Gaussian distribution was the best match. Concerning most of the northern sector, where the highest severity values were observed, the variable was best modeled by exponential distributions. Conversely, the central and intermediate regions of the NRB were mostly better adjusted by lognormal distributions.
The time series of D showed spatial patterns (Figure 5b) that resemble the average and maximum

Marginal Distributions
The multivariate distributions obtained by copula are generated from the marginal distributions of the individual variables. The advantage of the copula method is the flexibility to consider the probabilistic nature of each variable individually. Therefore, the marginal distributions are the probability distributions that best fit the variables under analysis. We tested the following probability distributions: normal (Gaussian), lognormal, exponential, logistic, Weibull and gamma [73]. According to the criterion of lowest AIC, the distributions that best fit the data are illustrated in Figure 5. The distributions selected for the severity series were the normal, lognormal and exponential ones (Figure 5a). The spatial pattern presented is similar to that observed in the fields of mean S and maximum S. Regarding the southeastern portion of the NRB, which had showed the lowest severity values (Figures 3a and 4a), the Gaussian distribution was the best match. Concerning most of the northern sector, where the highest severity values were observed, the variable was best modeled by exponential distributions. Conversely, the central and intermediate regions of the NRB were mostly better adjusted by lognormal distributions.
The time series of D showed spatial patterns (Figure 5b) that resemble the average and maximum D maps (Figure 4). Regarding the eastern portion, it was noted that the seasons with the longest drought durations were best modeled by the normal distribution.
Concerning the rest of the NRB, the duration of the drought was modeled by the exponential distribution.

Bivariate Joint Probability Distribution
After obtaining the marginal distributions for each NRB point, it is possible to choose which copula best fits the series under analysis. The chosen copula should be able to model sets of relations of a linear and nonlinear nature between the Severity and Duration series. Figure 6 shows the best joint distribution functions for each surface weather station chosen according to the lowest AIC value criterion. Throughout the NRB, two joint distributions were chosen, the Clayton and Frank copulas. The southeastern portion of the NRB has a cluster of modeled stations by copula of the Clayton-type [74]. This region, according to previous maps, is characterized by less severe and longer lasting droughts (Figures 3 and 4). Conversely, the northeastern sector of the NRB, characterized by higher severity and duration values (Figures 3 and 4), exhibits a cluster of Frank-type copula-shaped stations. The western side of the NRB exhibits no significant pattern on any dominant copula.

Bivariate Joint Probability Distribution
After obtaining the marginal distributions for each NRB point, it is possible to choose which copula best fits the series under analysis. The chosen copula should be able to model sets of relations of a linear and nonlinear nature between the Severity and Duration series. Figure 6 shows the best joint distribution functions for each surface weather station chosen according to the lowest AIC value criterion. Throughout the NRB, two joint distributions were chosen, the Clayton and Frank copulas. The southeastern portion of the NRB has a cluster of modeled stations by copula of the Clayton-type [74]. This region, according to previous maps, is characterized by less severe and longer lasting droughts (Figures 3 and 4). Conversely, the northeastern sector of the NRB, characterized by higher severity and duration values (Figures 3 and 4), exhibits a cluster of Frank-type copula-shaped stations. The western side of the NRB exhibits no significant pattern on any dominant copula.

Return Period
According to the selected joint distributions (Figure 6), it is possible to calculate the drought return period for a given level characterized by an S value (for Severity) and a D value (for Duration). Table 2 presents three drought classifications according to the definition of Wang et al, [33]. These classifications were used to illustrate the dynamics of the return period of the NRB seasons considering the hypothetical drought scenarios presented in Table 2.

Return Period
According to the selected joint distributions (Figure 6), it is possible to calculate the drought return period for a given level characterized by an S value (for Severity) and a D value (for Duration). Table 2 presents three drought classifications according to the definition of Wang et.al, [33]. These classifications were used to illustrate the dynamics of the return period of the NRB seasons considering the hypothetical drought scenarios presented in Table 2.
Applying the copula functions to the definitions of a weak drought (S = 3, D = 3), it is observed that the return period is longer in the western side of the NRB (Figure 7a), with an average interval of 17 months between one drought and another. Therefore, weak droughts are more frequent on the eastern side, as they have a shorter return period (~14 months). Regarding moderate droughts (S = 6, D = 6), there is some homogeneity over the NRB domain (Figure 7b). The western sector continues to show the longest return period (~26 months), but the difference is of only two months compared to the eastern sector (~24 months). Concerning severe droughts (S = 9, D = 9; Figure 7c), there is a clear inversion of the spatial pattern observed in previous maps. This time, the eastern portion has the longest return periods (~60 months) compared to the rest of the NRB domain (~45 months). This indicates that severe droughts are more frequent in the semiarid region. Table 2. Drought level parameters according to the criteria of Wang et.al, [33].

Drought Level Setting
Weak drought S = 3 D = 3 Moderate drought S = 6 D = 6 Strong drought S = 9 D = 9 Applying the copula functions to the definitions of a weak drought (S = 3, D = 3), it is observed that the return period is longer in the western side of the NRB (Figure 7a), with an average interval of 17 months between one drought and another. Therefore, weak droughts are more frequent on the eastern side, as they have a shorter return period (~14 months). Regarding moderate droughts (S = 6, D = 6), there is some homogeneity over the NRB domain (Figure 7b). The western sector continues to show the longest return period (~26 months), but the difference is of only two months compared to the eastern sector (~24 months). Concerning severe droughts (S = 9, D = 9; Figure 7c), there is a clear inversion of the spatial pattern observed in previous maps. This time, the eastern portion has the longest return periods (~60 months) compared to the rest of the NRB domain (~45 months). This indicates that severe droughts are more frequent in the semiarid region. Table 2. Drought level parameters according to the criteria of Wang et al, [33].

Drought Level Setting
The definitions of drought by [33] allow us to analyze drought dynamics from certain scenarios not necessarily representing actual values of droughts in the NRB. That is why we applied the average and maximum values of S and D obtained in the frequency analysis. This allowed us to obtain the empirical return period of the average droughts of each season and the return period of the most severe drought ever recorded in the time span comprised by the precipitation series. The definitions of drought by [33] allow us to analyze drought dynamics from certain scenarios not necessarily representing actual values of droughts in the NRB. That is why we applied the average and maximum values of S and D obtained in the frequency analysis. This allowed us to obtain the empirical return period of the average droughts of each season and the return period of the most severe drought ever recorded in the time span comprised by the precipitation series.
Regarding the average droughts observed in the NRB, the western portion showed the longest return period of approximately 26 months (Figure 8a), with peaks above 28 months spread over that region and in the north of the NRB. Some areas located in the center and southeast of the NRB exhibited minimum return periods, of about 20 months, whereas some locations of the states of AL, SE, RN, PB, CE and MA had similar return periods to that observed in the western sector. Regarding the most severe droughts (Figure 8b), the southern NRB stands out for its 90-150-month peaks. Added to the weather phenomena mentioned above that modulate climate variability in the eastern and northern NRB, there is also a significant role played by frontal systems on the coast of the state of Bahia [38]. The peak observed in that region indicates that the occurrence of droughts is related to Regarding the average droughts observed in the NRB, the western portion showed the longest return period of approximately 26 months (Figure 8a), with peaks above 28 months spread over that region and in the north of the NRB. Some areas located in the center and southeast of the NRB exhibited minimum return periods, of about 20 months, whereas some locations of the states of AL, SE, RN, PB, CE and MA had similar return periods to that observed in the western sector. Regarding the most severe droughts (Figure 8b), the southern NRB stands out for its 90-150-month peaks. Added to the weather phenomena mentioned above that modulate climate variability in the eastern and northern NRB, there is also a significant role played by frontal systems on the coast of the state of Bahia [38]. The peak observed in that region indicates that the occurrence of droughts is related to less incursions of frontal systems and a lower performance of EWD. The northern portion of the NRB showed the shortest return periods (~45 months) and is, therefore, the region with the most frequent severe droughts. Concerning the southeastern and northeastern sectors, the return periods are around 55 months.
According to Marengo et al, [1], there are widespread records of droughts in the NRB since the 15th century. The droughts that occurred within the period analyzed in this paper were as follows: 1966, 1970, 1976, 1979-1981, 1982-1983, 1986-1987, 1992-1993, 1997-1998, 2001-2002, 2005, 2010 and 2012-2016 Based on that record, the average return period between each event was 43 months, which is consistent with the numbers found in this study (Figure 7).
Part of such drought variability in the NRB is often associated with the occurrence of low frequency phenomena, such as the El Niño-Southern Oscillation (ENSO) [1,75]. During the occurrence of the ENSO positive/negative phase, there is a tendency for decreasing/increasing rainfall accumulations in the NRB. Another important low-frequency mode that modulates the NRB drought variability is the Tropical Atlantic Surface Temperature. The latitudinal annual migration of the ITCZ varies as a function of the surface temperature gradient between the Tropical North Atlantic and Tropical South Atlantic, referred to in literature as the Atlantic Dipole (AD) [17,76,77]. Namias et al, [15] further proposed that there are teleconnections between ENSO and AD that occur across the mid and high latitudes of the northern hemisphere through the North Atlantic Oscillation (NAO) patterns, as well as other modes of variability, such as the Pacific Decadal Oscillation (PDO) and the Atlantic Multidecadal Oscillation (AMO), described in detail by Kayano et al, [78,79].
Future scenario analysis for droughts in the NRB shows increased surface temperature and decreased precipitation [1]. The results of this study also imply an increase in the duration and frequency of droughts. Vieira et al, [80] pointed to an increase of approximately 5% of areas susceptible to desertification in the NRB from 2000 to 2010. Franchito et al, [81] indicate a 20% increase in aridity in the semiarid NRB region under climate change scenarios over the next 100 years. Significant positive drought trends during some NRB rainy seasons from 1980 to 2015 were shown by da Rocha Júnior et al, [35].
The NRB drought return period calculated for each season is useful information for planning long-term strategies to mitigate its effects [33,66]. The adjusted copula functions can be used operationally to monitor drought events, quantify their dangerousness during their evolution, and make it possible to adopt strategic actions in real time [82].
Geosciences 2019, 9, x FOR PEER REVIEW 12 of 17 less incursions of frontal systems and a lower performance of EWD. The northern portion of the NRB showed the shortest return periods (~45 months) and is, therefore, the region with the most frequent severe droughts. Concerning the southeastern and northeastern sectors, the return periods are around 55 months. According to Marengo et.al, [1], there are widespread records of droughts in the NRB since the 15th century. The droughts that occurred within the period analyzed in this paper were as follows: 1966, 1970, 1976, 1979−1981, 1982−1983, 1986−1987, 1992−1993, 1997−1998, 2001−2002, 2005, 2010 and 2012−2016 Based on that record, the average return period between each event was 43 months, which is consistent with the numbers found in this study (Figure 7).
Part of such drought variability in the NRB is often associated with the occurrence of low frequency phenomena, such as the El Niño-Southern Oscillation (ENSO) [1,75]. During the occurrence of the ENSO positive/negative phase, there is a tendency for decreasing/increasing rainfall accumulations in the NRB. Another important low-frequency mode that modulates the NRB drought variability is the Tropical Atlantic Surface Temperature. The latitudinal annual migration of the ITCZ varies as a function of the surface temperature gradient between the Tropical North Atlantic and Tropical South Atlantic, referred to in literature as the Atlantic Dipole (AD) [17,76,77]. Namias et.al, [15] further proposed that there are teleconnections between ENSO and AD that occur across the mid and high latitudes of the northern hemisphere through the North Atlantic Oscillation (NAO) patterns, as well as other modes of variability, such as the Pacific Decadal Oscillation (PDO) and the Atlantic Multidecadal Oscillation (AMO), described in detail by Kayano et.al, [78,79].
Future scenario analysis for droughts in the NRB shows increased surface temperature and decreased precipitation [1]. The results of this study also imply an increase in the duration and frequency of droughts. Vieira et.al, [80] pointed to an increase of approximately 5% of areas susceptible to desertification in the NRB from 2000 to 2010. Franchito et.al, [81] indicate a 20% increase in aridity in the semiarid NRB region under climate change scenarios over the next 100 years. Significant positive drought trends during some NRB rainy seasons from 1980 to 2015 were shown by da Rocha Júnior et.al, [35].
The NRB drought return period calculated for each season is useful information for planning long-term strategies to mitigate its effects [33,66]. The adjusted copula functions can be used operationally to monitor drought events, quantify their dangerousness during their evolution, and make it possible to adopt strategic actions in real time [82].

Conclusions
The Brazilian Northeast is the most populous dry region in the world, with a marked climate variability that exposes the population to the risk of drought having strong environmental and social impacts. The need for drought studies using advanced research methods, especially in the context of climate change, is critical to improve crisis management demanded by long periods of drought.
The most severe and frequent droughts observed in the NRB, on average, are concentrated in the northern portion. Seen in this region, the observed duration ranged from three and a half to five and a half months. Conversely, the sector that experienced the most severe droughts observed during the period 1961−2017 was the eastern NRB. This area, which extends along the coast of six states from RN to BA, had droughts which lasted from 14 to 24 months.
The marginal distributions that best modeled the S and D series were exponential, normal, and lognormal. Regarding the S series, it was observed that in the north there was the predominance of a better fit for exponential distributions, while in the central and eastern NRB, the lognormal and normal distributions were found to be better suited to deal with the data. Concerning the D series, a better performance of the lognormal distribution was noticed for nearly all of the NRB, except in the eastern sector where there was a higher frequency of normal distribution. Considering the choice of the marginal distributions of the D and S series for each season, the Frank and Clayton copula functions were chosen to model the drought return period in the NRB.
The analysis of scenarios with parameters of S and D given by [33] shows that, regarding mild droughts in the western portion of the NRB, they have longer return periods (~17 months) relative to the eastern side (~12 months). Regarding moderate droughts, there is a weak gradient pointing to the southwestern portion of the NRB. Moderate droughts throughout the NRB domain on average have a return period of around 22 months. Strong droughts, in their turn, presented opposite patterns to the ones found in weak droughts. The eastern side had a longer return period for severe droughts (~55 months) compared to the western area (~43 months).
When analyzing the return period based on observed S and D values, it was noted that, for average values, the return period in the entire region is around 24 months. Still in this scenario of average S and D, the southeast of the NRB stands out for presenting a return period of only about 20 months, whereas pockets above 28 months were found in the western and northern regions. Regarding maximum observed values of D and S, the eastern portion stands out for having the longest return period (>60 months), especially in the southeastern region, where a return period above 90 months was observed. The northern NRB had the shortest return period (~45 months), indicating that it is the NRB region with the highest frequency of severe droughts.
The results presented in this work are added to a hall of information useful to decision makers in terms of mitigating the negative effects of drought in the NRB. Moreover, they provide an initial methodological basis for an operational drought monitoring tool, contributing to real-time strategic decisions. Also, such knowledge can be applied in the quantification of the severity of future droughts from the representation of the hydrological cycle of models of the terrestrial system, as is the case of the Brazilian Earth System Model currently under development.