Copula-Based Multivariate Frequency Analysis of the 2012–2018 Drought in Northeast Brazil

The 2012–2018 drought was such an extreme event in the drought-prone area of Northeast Brazil that it triggered a discussion about proactive drought management. This paper aims at understanding the causes and consequences of this event and analyzes its frequency. A consecutive sequence of sea surface temperature anomalies in the Pacific and Atlantic Oceans, at both the decadal and interannual scales, led to this severe and persistent drought. Drought duration and severity were analyzed using run theory at the hydrographic region scale as decision-makers understand impact analysis better at this scale. Copula functions were used to properly model drought joint characteristics as they presented different marginal distributions and an asymmetric behavior. The 2012–2018 drought in Ceará State had the highest mean bivariate return period ever recorded, estimated at 240 years. Considering drought duration and severity simultaneously at the level of the hydrographic regions improves risk assessment. This result advances our understanding of exceptional events. In this sense, the present work proposes the use of this analysis as a tool for proactive drought planning.

Droughts have been reported in NEB since the colonial period [3,11,12]. Historically, the States of Ceará, Rio Grande do Norte, Paraiba, and Pernambuco are concentrated drought hotspots [13]. Drought hazards have caused massive migration and significant population death, such as the drought of 1877-1879, with human drought-related deaths estimated to be around 500,000 persons in Ceará State alone [12]. The population had no warning alert, and countless citizens chose to endure with minimal provisions before migrating to less impacted areas. Many perished in the process.
The inherent characteristics of the semi-arid region (strong seasonality coupled with high rainfall and discharge variability), shallow soils (most above the crystalline rock basement), and elevated evapotranspiration rates amplified drought-related impacts in NEB.
Societies adapt to the environment shaped by climate factors [12] and, in the case of Ceará, this adaptation was centered on two pillars. The first was the construction of large dams, used as public in the area. Section 3 investigates the large-scale dynamics that caused the 2012-2018 drought event and the main consequences derived from this event. Section 4 presents the data and methods used to define and model drought and its characteristics. Section 5 introduces the results of these analyses, including the duration and severity of the drought initiated in 2012 for each hydrographic region, and their univariate and bivariate return period. Section 6 places a summary with concluding remarks and discussions of the main results in comparison with the results from other authors.

Study Area
The State of Ceará is located in the northern portion of the northeast of Brazil (NNEB) (Figure 1). The hydrographic regions' code numbers are defined as a function of latitude. The hydrographic regions located in the North are classified as HR01 to HR07, those located in the Central area are HR08 to HR10, and the Southern hydrographic regions are HR11 and HR12. Ceará has more than 90% of its territory located in the Brazilian semiarid region (characterized by low precipitation levels, less than 800 mm per year, high evaporation rates, and shallow soils). Such characteristics make the region remarkably vulnerable to droughts. This vulnerability is increasing due to permanent land degradation, which puts 94% of NEB into moderate to high risk of desertification [7].
Water 2019, 11, x FOR PEER REVIEW 3 of 23 mechanisms causing rainfall in the area. Section 3 investigates the large-scale dynamics that caused the 2012-2018 drought event and the main consequences derived from this event. Section 4 presents the data and methods used to define and model drought and its characteristics. Section 5 introduces the results of these analyses, including the duration and severity of the drought initiated in 2012 for each hydrographic region, and their univariate and bivariate return period. Section 6 places a summary with concluding remarks and discussions of the main results in comparison with the results from other authors.

Study Area
The State of Ceará is located in the northern portion of the northeast of Brazil (NNEB) ( Figure  1). The hydrographic regions' code numbers are defined as a function of latitude. The hydrographic regions located in the North are classified as HR01 to HR07, those located in the Central area are HR08 to HR10, and the Southern hydrographic regions are HR11 and HR12. Ceará has more than 90% of its territory located in the Brazilian semiarid region (characterized by low precipitation levels, less than 800 mm per year, high evaporation rates, and shallow soils). Such characteristics make the region remarkably vulnerable to droughts. This vulnerability is increasing due to permanent land degradation, which puts 94% of NEB into moderate to high risk of desertification [7]. The rainy season in Ceará State is characterized by a distinct seasonality, extending from December to July. The main rainy season occurs from February to May, depending on oceanic and atmospheric conditions, and the principal mechanism that influences rainfall in this period is the Intertropical Convergence Zone (ITCZ) [12,14]. When the difference of the sea surface temperature (SST) of tropical north and south Atlantic, i.e., the Interhemispheric Tropical Atlantic Gradient The rainy season in Ceará State is characterized by a distinct seasonality, extending from December to July. The main rainy season occurs from February to May, depending on oceanic and atmospheric conditions, and the principal mechanism that influences rainfall in this period is the Intertropical Convergence Zone (ITCZ) [12,14]. When the difference of the sea surface temperature (SST) of tropical north and south Atlantic, i.e., the Interhemispheric Tropical Atlantic Gradient (IHTAG), is weaker, the ITCZ reaches its southernmost position, which usually occurs around March-April [14,44,45]. The interannual climatic variability of the NNEB is highly modulated by thermodynamic patterns that occur over the tropical Pacific and Atlantic Oceans. El Niño South Oscillation (ENSO) and IHTAG can perturb Walker and Hadley cells, causing drifts and consequently changing the intensity and period of the rainy season in the area [14,15,[45][46][47].
Depending on the intensity and period of the year in which it occurs, the ENSO warm phase is responsible for displacing the descending portion of the Walker cell. This phenomenon causes a zone of high pressure over NEB, which makes cloud formation in the region difficult and, consequently, influencing in years considered dry or very dry in the region [15,44,45]. The IHTAG, on the other hand, is capable of causing drought in the region when abnormal warming of tropical North Atlantic SST occurs, which creates a low-pressure zone in this part of the Atlantic Ocean and attracts the ITCZ towards North, avoiding precipitation over South American continent [45,46,48,49]. In addition, the association of the positive phase of North Atlantic SST concomitant with an El Niño event provides accentuated regional impacts on the climatic condition [50].
Despite the high-frequency interannual variability, there is also decadal climatic variability that can be influenced by low-frequency modes of SST anomalies. The pacific SST varies at a decadal time scale, a mode of SST variability known as Pacific Decadal Oscillation (PDO).
Kayano and Andreoli [47] found a significant climate teleconnection between the precipitation in NNEB and PDO. The Atlantic Ocean also has its low-frequency variability mode, referred to as Atlantic Multidecadal Oscillation (AMO) [51,52]. Linkages between AMO and seasonal climate variability over NNEB have also been found [53][54][55][56][57]. Periods with simultaneously positive (or negative) PDO and AMO phases result in a more predictable behavior of rainfall in the region, with values below (or above) than normally expected [58].

The 2012-2018 Drought
The shock of the 2012-2018 drought raised awareness in Ceará's management community to the importance of proactive drought management, primarily due to the severity of its impacts. It was also an opportunity to better comprehend the complex interactions between ocean and atmosphere and the consequences in the rainfall over NNEB. The leading causes of this event were related to the serial combination of high and low-frequency anomalies in SSTs that caused its persistence. Therefore, the 2012-2018 event can be understood as a series of consecutive one-year droughts.
Although El Niño is usually associated with dry periods and La Niña with wet periods in the NNEB, this drought started under the influence of a La Niña event. Rodrigues and McPhaden [59], analyzing how the 2011-2012 La Niña event could have caused the drought in 2012, found two different types of La Niña event: (1) the cooling concentrated in the eastern Pacific, causing a cooling of the tropical North Atlantic and warming of the tropical South Atlantic; and (2) the cooling concentrated in the central Pacific, causing the opposite SST gradient in the IHTAG. The first type, the classical understanding of La Niña, can bring rain to the NNEB. The second is the one that caused the drought of 2012, which induced migration of ITCZ towards the north. ENSO is a complex phenomenon and the full comprehension of its interactions with precipitation over the studied area is still in progress as new events occur. Additionally, an upper-level convergence over NEB associated with an upper-level divergence in Amazonia during the main rainy season in 2012 contributed to this drought [60].
The years of 2013 and 2014 did not present explicit forcing in the inter/intra annual scales, and the AMO/PDO process may have influenced the drought in these years. In 2013, the ENSO phenomenon presented conditions of neutrality; although, ITCZ operated north of its climatological position in response to the near-neutral but still warming condition of surface waters in the tropical North Atlantic. The ITZC position, combined with westward anomalous displacement of humidity at high levels, contributed to rainfall below average in NNEB.
In 2014, neither El Niño nor IHTAG presented strong signals, and spatial variability of rainfall anomalies was found in NEB. However, a climatic condition that may have contributed to drought during the period was an anticyclonic anomaly detected in southeastern Brazil and considered one of the most critical factors of the 2014-2015 drought that affected Southeast Brazil [61]. This system had an extension to NEB, affecting the area since 2012 [2].
In 2015, the expansion of positive SST anomalies along the equatorial Pacific Ocean indicated the full establishment of the ENSO phenomenon. This El Niño persisted and gained force during 2015 and 2016, influencing the below-average rainfall in 2016 and 2017. In 2017, the El Niño condition retreated, and a La Niña configuration initiated. This state was favorable to indicate the end of the drought, but the warming of the North Tropical Atlantic Ocean was also detected, negatively influencing the rain in NEB.
In 2018, IHTAG indicated a negative phase, especially around the end of the rainy season, and, in association with La Niña configuration over the equatorial Pacific, contributed to rainfall around the climatological average over NNEB. This configuration provided enough rainfall to recover from the drought state in the majority of NNEB; however, few places still present persistence of this event. Figure 2 presents the time series of the cumulative rainfall anomaly and main climatic indexes that present teleconnections to precipitation in the region, NINO 3.4, IHTAG, PDO and AMO for the period of 2009 to 2018 (https://www.esrl.noaa.gov/psd/data/climateindices/list/). The accumulated rainfall deficit of the 2012-2018 drought is 1225 mm, 1.5 times the yearly climatological rainfall, 800.6 mm.
Water 2019, 11, x FOR PEER REVIEW 5 of 23 during the period was an anticyclonic anomaly detected in southeastern Brazil and considered one of the most critical factors of the 2014-2015 drought that affected Southeast Brazil [61]. This system had an extension to NEB, affecting the area since 2012 [2]. In 2015, the expansion of positive SST anomalies along the equatorial Pacific Ocean indicated the full establishment of the ENSO phenomenon. This El Niño persisted and gained force during 2015 and 2016, influencing the below-average rainfall in 2016 and 2017. In 2017, the El Niño condition retreated, and a La Niña configuration initiated. This state was favorable to indicate the end of the drought, but the warming of the North Tropical Atlantic Ocean was also detected, negatively influencing the rain in NEB.
In 2018, IHTAG indicated a negative phase, especially around the end of the rainy season, and, in association with La Niña configuration over the equatorial Pacific, contributed to rainfall around the climatological average over NNEB. This configuration provided enough rainfall to recover from the drought state in the majority of NNEB; however, few places still present persistence of this event. Figure 2 presents the time series of the cumulative rainfall anomaly and main climatic indexes that present teleconnections to precipitation in the region, NINO 3.4, IHTAG, PDO and AMO for the period of 2009 to 2018 (https://www.esrl.noaa.gov/psd/data/climateindices/list/). The accumulated rainfall deficit of the 2012-2018 drought is 1225 mm, 1.5 times the yearly climatological rainfall, 800.6 mm.  The long duration and severity of the 2012 drought caused many impacts on the socio-economy and environment in NEB. Figure 3 presents the stored water volume per hydrographic region and the total accumulation in the state of Ceará. The total water storage capacity of Ceará State is 18,500 hm 3 (26% in the north, 57% in central and 17% in the south region). This storage water decreased by 63%, from 131 × 10 6 hm 3 in December 2011 to 13 × 10 6 hm 3 by the end of 2016, with some hydrographic regions with total collapse. This issue was more accentuated and prolonged in central and southern regions, i.e., HR8-HR12, which represent 74% of the state's total accumulation capacity. With the prolongation of the drought, the small and medium reservoirs (with storage capacity below 75 hm 3 according to Ceará State's declaration no. 23.068/1994 [62]) started to collapse, both in terms of quantity and quality, enhancing the costs of capturing and distributing water at longer distances. In red are the accumulated rainfall deficit and oceanic index associated with dry years in the region. In blue, the exceeding rainfall and the oceanic index associated with wet years.
The long duration and severity of the 2012 drought caused many impacts on the socio-economy and environment in NEB. Figure 3 presents the stored water volume per hydrographic region and the total accumulation in the state of Ceará. The total water storage capacity of Ceará State is 18,500 hm³ (26% in the north, 57% in central and 17% in the south region). This storage water decreased by 63%, from 131 × 10 6 hm³ in December 2011 to 13 × 10 6 hm³ by the end of 2016, with some hydrographic regions with total collapse. This issue was more accentuated and prolonged in central and southern regions, i.e., HR8-HR12, which represent 74% of the state's total accumulation capacity. With the prolongation of the drought, the small and medium reservoirs (with storage capacity below 75 hm³ according to Ceará State's declaration no. 23.068/1994 [62]) started to collapse, both in terms of quantity and quality, enhancing the costs of capturing and distributing water at longer distances. The water shortage also affected the water quality of those reservoirs, especially regarding eutrophication and an increase in the concentration of salts due to the low inflow periods, higher evaporation, and anthropogenic actions. Santos et al. [63], for instance, monitored the water quality of the biggest reservoir in Ceará State (i.e., the Castanhão, located at HR10) from November 2011 to May 2014 and found that it went from an initial oligotrophic condition, i.e., low nutrient values, in 2011 to an eutrophication condition, i.e., high nutrient values, with the decrease in its accumulated volume by 2014. Increased water treatment costs arose from this condition.
In response to this water shortage, federal and state actions were taken to build a series of emergency pipelines, drill wells, construct water cisterns and distribute water through water tank trucks to meet the demands in rural and urban areas in Ceará. The reactive characteristic of the measures taken is implied in increased associated costs, as no previous planning for these actions existed.
The impacts on the state's agriculture were felt at different timescales, depending on the type of agriculture used. In the first two years of the drought, 2012-2013, rainfed agriculture was strongly impacted, and many farmers completely abandoned their cultures. The abandoned soil enabled natural vegetation, adapted to dry conditions, to recover, even during prolonged drought periods. Irrigated agriculture, on the other hand, suffered practically no impact at the start of the drought, since the large multi-annual reservoirs guaranteed its supply. Those reservoirs initiated the drought  The water shortage also affected the water quality of those reservoirs, especially regarding eutrophication and an increase in the concentration of salts due to the low inflow periods, higher evaporation, and anthropogenic actions. Santos et al. [63], for instance, monitored the water quality of the biggest reservoir in Ceará State (i.e., the Castanhão, located at HR10) from November 2011 to May 2014 and found that it went from an initial oligotrophic condition, i.e., low nutrient values, in 2011 to an eutrophication condition, i.e., high nutrient values, with the decrease in its accumulated volume by 2014. Increased water treatment costs arose from this condition.
In response to this water shortage, federal and state actions were taken to build a series of emergency pipelines, drill wells, construct water cisterns and distribute water through water tank trucks to meet the demands in rural and urban areas in Ceará. The reactive characteristic of the measures taken is implied in increased associated costs, as no previous planning for these actions existed.
The impacts on the state's agriculture were felt at different timescales, depending on the type of agriculture used. In the first two years of the drought, 2012-2013, rainfed agriculture was strongly impacted, and many farmers completely abandoned their cultures. The abandoned soil enabled natural vegetation, adapted to dry conditions, to recover, even during prolonged drought periods. Irrigated agriculture, on the other hand, suffered practically no impact at the start of the drought, since the large multi-annual reservoirs guaranteed its supply. Those reservoirs initiated the drought with elevated accumulated levels regarding the previous rainy year of 2011. With the persistence of drought and the consequent decrease in accumulated levels, the reduction and posterior interruption of water use permits for irrigation were determined to save water for the prioritized human water supply, according to Brazilian water law.
In this sense, a series of crises management measures were promoted by Federal programs, such as: Programa Garantia Safra, granted to farmers that lost at least 50 percent of their production; Bolsa Estiagem, that distributed US$40/month for smallholder farmers; subsidized prices for selling maize to feed animals; expansion of emergency credit lines for farmers, traders and industry sectors; and revised debt of farmers [1]. Despite the devastating impacts on agricultural, livestock, and industrial activities, this extreme drought did not lead to human losses nor migration as happened in the past; this lower social disturbance is associated with government social programs [1,12].

Data
The series of daily precipitation from 1911 to 2018 used to analyze drought in Ceará State was obtained from the Brazilian National Water Agency (http://www.snirh.gov.br/hidroweb/). The average areal rainfall for each of the twelve hydrographic regions of the Ceará State was obtained by interpolating the daily precipitation at each rain gauge according to the inverse distance weighting (IDW) method with exponent two into grid points with 0.05 • size. Further, the average of the interpolated values was extracted for those inside each analyzed area. The use of hydrographic region' scale tends to reduce random fluctuations of a point-based approach, homogenizing drought analysis, and reinforcing correlation with socio-economic impacts felt in this planning scale.

Drought Analysis
The drought analysis was based on the calculation of the Standard Precipitation Index (SPI) [64] with 12 months aggregated timescale (SPI12). The time scale used for the calculation of the index is directly related to the time required for the effects of drought to be felt on the different activity sectors and the region's water resources [65]. A new time-series was created with only December SPI12 values (SPI12 DEC ) to represent the accumulated annual information. This discretization process was performed to remove the continuous information provided by SPI moving window. By doing so, the objective was to archive independent random variables that represent the total annual precipitation, smoothing the temporal series and avoiding spurious information of SPI influenced by above or below-average precipitation in months of the dry season (July to December).
The calculation of drought duration and severity characteristics for each hydrographic region was obtained through run theory, as proposed by Yevjevich [66]. Each drought event is defined as the proportion of time all values of a variable X t are below a selected truncation level. Specifically for SPI, the duration of a given drought event is determined by summing the periods that this event remained below a certain threshold, in this paper, X t = 0. Shiau [24] used this threshold as it avoids the division of spurious droughts that occurs inside one longer drought, which makes sense as social and environmental impacts are more significant in prolonged droughts than in consecutive shorter droughts. Figure 4 illustrates this process. Drought events 1, 2, and 3 are orange. The severity is given by the summation of SPI values during one event, according to Equation (1): where S is the severity, and D is the duration, and SPI12 DEC is the SPI value discretized for every December considering the aggregated time-period of 12 months.

Statistical Inference
Once the drought duration and severity characteristics were separated from the original time series by run theory application, data analysis could be performed to deduce the properties of the variables samples and adjust to a population. The distribution function that better represents drought duration and severity has not yet established an agreement. Exponential and Gamma were proposed by Shiau [24] for univariate modeling the drought characteristics; however, many authors prefer to perform a goodness-of-fit test to find the families that best represent the analyzed sample for each region and drought [31,67].
Thus, both duration and severity time series were adjusted for the univariate Log-normal, Exponential, Weibull, Gamma, and Logistic probability distribution families. The parameters were chosen based on the maximum likelihood estimation (MLE) method. The Akaike information criteria (AIC) indicated the candidate distributions that best fitted the data. AIC is a parsimonious estimator of the relative quality of statistical models that penalizes overfitting.
Regarding the bivariate model, this study focused on the use of copula functions to model the dependence structure among marginal distribution functions. The bivariate joint distribution function ( , ) , where and are the random correlated variables, drought duration, and severity, with respective marginal distributions ( ) and ( ), is given by the copula function [ ( ), ( )], according to the Equation (2): where ( ) and ( ) are equal to and , with , ∈ (0,1). The copula functions can be classified as Meta-elliptic and Archimedean copulas: the first is symmetric, presenting no tail dependence; the second is more flexible and can present upper or lower tail dependence. In this study, three Archimedean copulas, Clayton, Frank, and Gumbel, and two Meta-elliptic copulas, Gaussian and t-Student, were used as candidates to identify the family that was best suited to model the dependence structure between the duration and severity. Equations (3)-(7) present the formulations of the candidate copula families, where , ρ, and v are the copula function parameters. Clayton Frank

Statistical Inference
Once the drought duration and severity characteristics were separated from the original time series by run theory application, data analysis could be performed to deduce the properties of the variables samples and adjust to a population. The distribution function that better represents drought duration and severity has not yet established an agreement. Exponential and Gamma were proposed by Shiau [24] for univariate modeling the drought characteristics; however, many authors prefer to perform a goodness-of-fit test to find the families that best represent the analyzed sample for each region and drought [31,67].
Thus, both duration and severity time series were adjusted for the univariate Log-normal, Exponential, Weibull, Gamma, and Logistic probability distribution families. The parameters were chosen based on the maximum likelihood estimation (MLE) method. The Akaike information criteria (AIC) indicated the candidate distributions that best fitted the data. AIC is a parsimonious estimator of the relative quality of statistical models that penalizes overfitting.
Regarding the bivariate model, this study focused on the use of copula functions to model the dependence structure among marginal distribution functions. The bivariate joint distribution function H(d, s), where D and S are the random correlated variables, drought duration, and severity, with respective marginal distributions F D (d) and F S (s), is given by the copula function C[F D (d), F S (s)], according to the Equation (2): where F D (d) and F S (s) are equal to u and v, with u, v ∈ (0, 1). The copula functions can be classified as Meta-elliptic and Archimedean copulas: the first is symmetric, presenting no tail dependence; the second is more flexible and can present upper or lower tail dependence. In this study, three Archimedean copulas, Clayton, Frank, and Gumbel, and two Meta-elliptic copulas, Gaussian and t-Student, were used as candidates to identify the family that was best suited to model the dependence structure between the duration and severity. Equations (3)-(7) present the formulations of the candidate copula families, where θ, ρ, and v are the copula function parameters. Frank Gumbel exp − (− ln u 1 ) θ + (− ln u 2 ) θ 1 θ Gaussian Water 2020, 12, 834 9 of 22 The inference function from margins (IFM) method [68] was used to estimate copula parameters. IFM is a parametric method that consists of the previous definition of marginal distributions used to transform samples in the (0,1) interval. Thus, transformed samples are jointly modeled by estimating the parameters for the families of the candidate copula using the maximum likelihood method. The minimum value of AIC was used to find the best fit model around the candidate copulas. Brechmann and Schepsmeier [69] defined the AIC relationship with a bivariate copula model and its respective parameter (θ), according to Equation (8).
where i = 1, . . . , N are the observations of the variables modeled and k the number of estimated parameters in the model.

Frequency Analysis
To better prepare for the occurrence of future droughts, one analysis that can be integrated into the risk management is the estimation of return periods of past drought events through a process known as frequency analysis. Independent and stationary time series are needed to perform the frequency analysis [20].

Univariate Return Period
The calculation of the return period represents the expected period between the occurrence of two events with the same or superior magnitude. The return period of drought duration (T D ) and severity (T S ) are described as a function of the expected interarrival time E(L) and the cumulative distribution functions (CDF) of the drought characteristic F D (d) and F S (s), as expressed in Equations (9) and (10) [24,67,70,71]. For the return period of a time series with annual recurrences, such as annual maxima precipitation, the E(L) is equal to one. However, droughts are not supposed to occur every year, and E(L) is found by estimating the mean value between the occurrences of droughts.

Bivariate Return Period
According to Shiau [24], the joint return period of duration and severity can be defined in two cases: the return period for D ≥ d or S ≥ s and return period for D ≥ d and S ≥ s. Both definitions of joint return period for copula-based drought events are described by Equations (11) and (12), respectively: where T DorS is the return period for D ≥ d or S ≥ s; T D&S is the return period for D ≥ d and S ≥ s.

Drought Analysis
Droughts areby definition, extreme events and any proactive measure must be previously defined based on magnitude and frequency of occurrence. In drought-prone areas such as the Brazilian semi-arid region, the high interannual variability of hydrologic conditions must be considered on the standard operational routine, and only exceptional events should justify special treatment and institutional intervention [72]. Therefore, a scientific criterion is required to quantify the frequency of each event. Figure 5 indicates the SPI calculated for the mean precipitation over the Ceará State for the aggregated time-period of 1 to 35 months. The darker colors in 2012-2013 indicate that this was the most critical period of the analyzed drought, and the following individual years were smoothed. Therefore, the gravity was the combination of its strong beginning with the abnormal sequence of dry years.
where is the return period for ≥ ≥ ; & is the return period for ≥ and ≥ .

Drought Analysis
Droughts areby definition, extreme events and any proactive measure must be previously defined based on magnitude and frequency of occurrence. In drought-prone areas such as the Brazilian semi-arid region, the high interannual variability of hydrologic conditions must be considered on the standard operational routine, and only exceptional events should justify special treatment and institutional intervention [72]. Therefore, a scientific criterion is required to quantify the frequency of each event. Figure 5 indicates the calculated for the mean precipitation over the Ceará State for the aggregated time-period of 1 to 35 months. The darker colors in 2012-2013 indicate that this was the most critical period of the analyzed drought, and the following individual years were smoothed. Therefore, the gravity was the combination of its strong beginning with the abnormal sequence of dry years.  Source: FUNCEME (2019). Drought events for the twelve hydrographic regions of Ceará State were identified using SPI12 DEC . For the sake of simplicity, SPI12 DEC will be treated by SPI from now on. The SPI values for the 2011-2018 drought exceeded the extreme condition (threshold equals −2.0, according to Mckee, 1993) for the majority of hydrographic regions during the drought onset in 2012. However, it was not the first time a drought with such magnitude had occurred in the region, such as the events around 1920, and during the 1950s decade, as shown in Figure 6. The columns in Figure 6 also show the spatial coverage of extreme events, such as the warm colors at the dry years of 1915,1919,1932,1958,1983,1993 and 2012, and the cold colors at the wet years of 1917, 1924, 1964, 1974 and 1985. This fact is associated with large-scale systems such as ENSO teleconnections and Atlantic circulation. This regional behavior is detrimental to drought management as all basins are uniformly affected, making it challenging to transfer water between hydrographic regions.
From Figure 6, it is also possible to see that the drought that started in 2012 has different ending times. For the hydrographic regions closer to the ocean, it ended between 2016 and 2017, indicated by light green colors. For those regions located more centrally and in the southern regions, the drought persisted until 2018, with no clear definition of ending for some of them yet. Also, consecutive drought years such as the analyzed one and covering almost all hydrographic regions can present enormous negative impacts and be detected in 1930-1933, 1941-1943, 1951-1956, 1979-1983 and 1990-1993, showing the high climatic variability existent in the region. From Figure 6, it is also possible to see that the drought that started in 2012 has different ending times. For the hydrographic regions closer to the ocean, it ended between 2016 and 2017, indicated by light green colors. For those regions located more centrally and in the southern regions, the drought persisted until 2018, with no clear definition of ending for some of them yet. Also, consecutive drought years such as the analyzed one and covering almost all hydrographic regions can present enormous negative impacts and be detected in 1930-1933, 1941-1943, 1951-1956, 1979-1983 and 1990-1993, showing the high climatic variability existent in the region. The analysis of the descriptive statistics of drought events (Table 1) showed that the hydrographic regions presented between 22 and 26 drought events over the period 1911-2018. In this period, a drought occurred every 4 to 5 years in each hydrographic region in Ceará State, as shown by the average inter-arrival time. In general, in the North, there were more droughts; however, they were shorter in duration and less severe than in the central and southern regions, where generally The analysis of the descriptive statistics of drought events (Table 1) showed that the hydrographic regions presented between 22 and 26 drought events over the period 1911-2018. In this period, a drought occurred every 4 to 5 years in each hydrographic region in Ceará State, as shown by the average inter-arrival time. In general, in the North, there were more droughts; however, they were shorter in duration and less severe than in the central and southern regions, where generally fewer droughts happened, but these were longer-lasting and more severe. The longest drought occurred in region 08, lasting 10 years, and the most severe in region 11, both located in the central and southern regions. HR08 presents the highest coefficient of variation (CV) for both duration and severity, indicating that this hydrographic region has the highest exposure to extreme droughts. Table 1 also shows that, for most hydrographic regions, the current drought is the most severe and prolonged ever recorded.
Although this analysis indicated that the 12 hydrographic regions of Ceará State have similar univariate descriptive statistics, the dependence structure between modeled variables dictates the joint behavior that is the object of analysis in this paper. To be able to model the joint distribution, the construction of the marginal distribution approach was used. Using the AIC as decision criteria to choose best-fit distributions, the duration series were best modeled by Log-normal distribution, while an Exponential distribution better represented the severity series.  Figure 7 shows the scatterplot of drought severity and duration. The 2012-2018 drought is one of the most adverse events ever recorded for most hydrographic regions, being compared to the droughts of 1951-1956 and 1978-1983. The dependence structure of drought duration and severity presented a tendency to become narrowly correlated with the increase of the values showing an upper tail correlation (Figure 7). Thus, a simple linear regression model hardly models this kind of asymmetric correlation. Copula functions, however, can meet this type of dependence structure. The drought initiated in 2012 is plotted in red. This analysis agrees with the one presented in Table 1, putting this event as one of the most severe ever recorded for all the hydrographic regions, which associated with population growth, increased water consumption and the reactive measures taken may explain the massive impacts caused by this drought. By knowing the marginal distributions, the copula functions could be fitted to the data. The Survival Clayton (180° rotated Clayton) and Gumbel, both asymmetric Archimedean copulas, were chosen. Table 2 shows the marginal distribution function, the copula functions, their respective By knowing the marginal distributions, the copula functions could be fitted to the data. The Survival Clayton (180 • rotated Clayton) and Gumbel, both asymmetric Archimedean copulas, were chosen. Table 2 shows the marginal distribution function, the copula functions, their respective parameters, and the Kendall's Tau correlation coefficient (τ) for each hydrographic region. The moderate τ values show that the founded duration and severity are not highly correlated. One possible reason for this is that the drought threshold equals zero, which selected a high number of droughts that lasted only one year. The asymmetric behavior of drought characteristics provided by this threshold is still able to be modeled by taking advantage of copula capabilities to model tail dependence. Also, this threshold is still adequate as it incorporates the impacts of drought enhancement and recovery. This moderate correlation additionally shows the importance of doing a multivariate analysis as drought duration does not necessarily indicate extreme severity.

Frequency Analysis
As stated by Haan [20], in order to perform frequency analysis, it is primarily necessary to check for independence and stationarity of the SPI time series. The independence was achieved by the discretization of the SPI12 into SPI12 DEC , as discussed in Section 4.2. Figure 8 shows the autocorrelation to the time series of SPI12 and SPI12 DEC for the HR04 as an example. The SPI12 time series presented strong serial dependence, i.e., values at some time t are statistically dependent to other lagged values, due to the moving window used to compute SPI12 values. The discretization performed by using SPI12 DEC time series provided the required independence to perform frequency analysis.
hypothesis of no trend was tested against the alternative hypothesis of the monotonic trend (not shown). From the 12 hydrographic regions tested, only two, HR05 and HR09, presented a statistically significant downward trend at a significance level of 0.05 in the 12 time series. Therefore, for those regions, the return period can be expected to be overestimated, which means that more frequent events can happen. However, those regions have smaller population density and do not contribute to the water transfer systems that provide water to coastal areas with higher population densities. for HR04. The 12 time series presented strong autocorrelation due to the moving window used to compute its values. In order to provide frequency analysis, independence of the time series was archived by performing a discretization by using 12 .
Once the considerations to perform frequency analysis are analyzed, the joint distribution function modeled based on the marginals and using copula functions can calculate the return period Figure 8. The differences between autocorrelation of time series of SPI12 and the discretized SPI12 DEC for HR04. The SPI12 time series presented strong autocorrelation due to the moving window used to compute its values. In order to provide frequency analysis, independence of the time series was archived by performing a discretization by using SPI12 DEC The Mann-Kendall (MK) test [73,74] was used to detect the trends in SPI12 DEC time series for the hydrographic regions, which is commonly used in hydrology and meteorology [75]. The null hypothesis of no trend was tested against the alternative hypothesis of the monotonic trend (not shown). From the 12 hydrographic regions tested, only two, HR05 and HR09, presented a statistically significant downward trend at a significance level of 0.05 in the SPI12 DEC time series. Therefore, for those regions, the return period can be expected to be overestimated, which means that more frequent events can happen. However, those regions have smaller population density and do not contribute to the water transfer systems that provide water to coastal areas with higher population densities.
Once the considerations to perform frequency analysis are analyzed, the joint distribution function modeled based on the marginals and using copula functions can calculate the return period as indicated by Equations (7) and (8). Figure 9 presents the return periods for the hydrographic regions for both the "and" and "or" cases. The return periods are presented in the form of contour lines. Different combinations of drought duration and severities can provide the same value for the return period. In the "or" case, contour lines do not cross the axes, while horizontal and vertical axes bound the "and" cases. It can be seen that "and" cases have higher return periods than the "or" ones, as the first analysis is more restrictive than the second. The information provided by Figure 9 can also specify the return period of a given event by providing its duration and severity. This functionality enables its use in proactive drought plans as the return period of any given drought can easily be found by providing the associated duration and severity. The 2012 event is highlighted in red. Most of the drought events that occurred in all hydrographic regions have a return period below the 64-year isoline. More extreme events, such as 1930-1933, 1941-1943, 1951-1956, 1979-1983 and 1990-1993, in addition to the 2012-2018 event, have a more extended return period. Table 3 summarizes the information from return period for univariate, duration and severity, and multivariate, "or" and "and" cases for these droughts and rank those return periods with the set of events recorded. The comparison between univariate return periods (T D and T S ) of the 2012-2018 drought for all analyzed regions presented no clear pattern of which one presents the highest values. For HR05, HR08, HR09, and HR11, hydrographic regions located in the south, central and southern regions, T S > T D , for the others T D > T S .
It is important to observe that compound events must satisfy following inequalities: T DorS < min( T D , T S ), i.e., the compound return period for the "or" case must be inferior to the minimum of univariate return period of those drought characteristics. As it is a more permissive event, only one of the two conditions must be satisfied. Also, T D&S > max( T D , T S ) implies that the compound return period for the "and" case must be superior to the maximum of the univariate return periods. As it is a more restrictive event, both conditions must be satisfied. Therefore, the joint return periods for the "and" cases are consistently higher than the univariate approach. These results indicate the necessity to consider the joint relationship between drought characteristics to real represents its recurrence as the correlation between drought characteristics are proportional to its damage potential. The rank also shown in Table 3 put the 2012-2018 drought as one of the three highest, most exceptional droughts that occurred between 1911-2018 for all hydrographic regions. Also, although more extended droughts had occurred in some hydrographic regions, the severity of this drought is highlighted, indicating the importance of multivariate analysis of drought events.
Drought risk for any region is a product of the region's exposure to a predefined event and the vulnerability of society to this event [76]. The return period can express the drought exposure as it incorporates the probability of the occurrence of an event. Therefore, it is interesting to analyze the exposure to drought hazard in the different hydrographic regions of an event with similar characteristics to the 2012-2018 drought. Thus, the return period of an event with average characteristics of the analyzed drought in the 12 regions (duration equals 6 years and severity equals 6) was calculated ( Figure 10). It shows a clear distinction between northern regions with central and southern areas. The south is the region which is more susceptible to severe and persistent droughts, such as the 2012-2018 event, and it is where the main reservoirs are located. The north is less affected by long drought as it is affected by intra-annual variability caused by oceanity conditions. Also, it is less dependent on ITCZ position as even a slight modification of its climatological position can still provide precipitation to the area. On the contrary, precipitation in the central and southern regions is more dependent on ITCZ, and consequently, to IHTAG modulation. Drought risk for any region is a product of the region's exposure to a predefined event and the vulnerability of society to this event [76]. The return period can express the drought exposure as it incorporates the probability of the occurrence of an event. Therefore, it is interesting to analyze the exposure to drought hazard in the different hydrographic regions of an event with similar characteristics to the 2012-2018 drought. Thus, the return period of an event with average characteristics of the analyzed drought in the 12 regions (duration equals 6 years and severity equals 6) was calculated ( Figure 10). It shows a clear distinction between northern regions with central and southern areas. The south is the region which is more susceptible to severe and persistent droughts, such as the 2012-2018 event, and it is where the main reservoirs are located. The north is less affected by long drought as it is affected by intra-annual variability caused by oceanity conditions. Also, it is less dependent on ITCZ position as even a slight modification of its climatological position can still provide precipitation to the area. On the contrary, precipitation in the central and southern regions is more dependent on ITCZ, and consequently, to IHTAG modulation. The lower the return period, the higher the exposure to drought hazard.
The presented drought frequency analysis indicates the recurrence of an event with magnitude equal to or greater than the one of the 2012-2018 drought for each hydrographic region in Ceará. It indicates that the joint return period is always higher than the univariate approach. This result indicates the necessity to consider the joint characteristics to understand the real exceptionality of extreme events. Another impressive result was that the northern areas are less susceptible to exceptional droughts, such as the analyzed one. As the main reservoir storage capacity is localized in The lower the return period, the higher the exposure to drought hazard.
The presented drought frequency analysis indicates the recurrence of an event with magnitude equal to or greater than the one of the 2012-2018 drought for each hydrographic region in Ceará. It indicates that the joint return period is always higher than the univariate approach. This result indicates the necessity to consider the joint characteristics to understand the real exceptionality of extreme events. Another impressive result was that the northern areas are less susceptible to exceptional droughts, such as the analyzed one. As the main reservoir storage capacity is localized in central and southern regions, this indicates that Ceará's water reserves are concentrated in the more vulnerable areas to the occurrence of prolonged and severe droughts.

Discussions and Conclusions
The northeast of Brazil (NEB) has experienced one of its worst droughts ever recorded, from 2012 to 2018. The leading causes were associated with anomalies in SST in the equatorial Pacific and Atlantic oceans caused by decadal and interannual variability modes. The serial combination and association of the climatic phenomenon (i.e., La Niña with the cooling occurring at central Pacific, the prevalence of tropical North Atlantic warming, AMO/PDO low-frequency modulations and El Niño) influenced the ITCZ and the Walker Circulation Cell to inhibit the occurrence of precipitation over NNEB. In Ceará State, the accumulated rainfall deficit of the 2012-2018 drought was 1225 mm, 1.5 times the yearly climatological rainfall.
NEB is known as a drought-prone region with considerable adaptive capacity, both in terms of increased water infrastructure and management. This resilience is based on learned experiences, acquired from its drought antecedents. This capacity has recently been questioned due to the magnitude of the current drought and the emergency measures taken to cope with it. Those measures helped to mitigate social damage that historically occurred in the most extreme droughts, such as human losses and massive migration [12]; however, they were taken under a "reactive" management paradigm, which could not handle some of the higher economic losses suffered by Ceará State.
We believe that proactive drought management can deal with some of the issues not addressed by its "reactive" predecessor. As stated by Gutiérrez et al. [1], this drought was the trigger needed to start a discussion regarding proactive drought management in NEB. Institutional relations between different public bodies and forums discussing the topic of drought have improved their performance by establishing monitoring processes and by incorporating active memory to elaborate proactive drought plans. Frequency analysis should be used in proactive drought plans to preserve the memory of past events. Therefore, proactive drought plans should use frequency analysis as it enables a scientific identification of drought recurrence, which can be used as a preparation tool for the mitigation of future droughts.
The univariate approach to calculate the drought return period has traditionally dominated the drought frequency analysis and it is a common practice in Brazil [2,5,6,10,72,77]. Martins et al. [10] estimated the return period of the 2012-2016 drought that occurred on the most significant water system in NEB, São Francisco River Basin, as 135 years using the univariate approach. However, the impact caused by a drought event may vary according to its duration and severity. Although these characteristics are correlated, their associated behavior can provide synergetic impacts that could be missed by a univariate approach.
The framework of bivariate frequency analyses can represent the exceptionality of drought events as the correlation between droughts characteristics are proportional to its damage potential, i.e., the negative impacts associated with a short but extremely severe drought may be stronger than another longer but less severe drought. Copula functions were useful to accurately model the dependence structure of drought characteristics as they presented different marginal distributions and due to observed upper tail dependence in the joint behavior. Thus, Gumbel and Survival Clayton asymmetric Archimedean copulas were chosen.
The hydrographic region was chosen as the planning scale in line with Brazilian water law that states it as the territorial management unit. There are significant benefits to use this scale, including better representation of socio-economic and environmental relationship existent between water supply and demand (e.g., precipitation, runoff, water reserves accumulated in reservoirs, associated demands for agricultural and urban uses). Also, it benefits from the capability to consolidate information that could otherwise randomly fluctuate in a point-based analysis.
The 2012-2018 drought in Ceará State had the highest mean bivariate return period ever recorded, presenting long persistence, substantial severity and spatial coverage. The mean joint return period, considering the "and" case, was 240 years (maximum of 499 years in HR05). The mean univariate return period of the 2012-2018 drought for the 12 hydrographic regions located in Ceará State was 169 years for the duration (maximum of 465 years in HR07), and 141 years for drought severity (maximum of 275 years in HR11). The bivariate analysis consistently presented higher values than the univariate, indicating the necessity to consider the joint behaviors to avoid underestimation of drought impacts. Similar characteristics to this drought were presented earlier in the 1951-1956 and 1978-1983 events for some regions, with a mean joint return period of 145 and 135 years, respectively.
The severity of this drought was influenced by the first two years, 2012 and 2013, added to the long final sequence; although the devastating impacts suffered from the current drought, having started with the most severe part of the drought, served as a critical warning. This opportunity increased the capability to mitigate drought effects in the area, but early warning and monitoring systems must be prepared to anticipate actions in future droughts that may not start with the same severity. Most of the other events presented a bivariate return period of less than 64 years for all hydrographic regions. Ceará State is more likely to present another drought with the same characteristics as the one here analyzed than California to a drought that occurred over the same period, 2012-2015, which has a return period estimated as 1400 years by Kwon and Lall [67]. This fact reflects the extreme variability and frequent drought recurrence existent in the region, showing the necessity to cope drought events with state-of-art techniques proper.
Knowing the exposure to drought is fundamental when planning measures to mitigate drought. The analysis presented here can inform decision-makers as to which areas are more susceptible to the occurrence of future droughts. This analysis indicated that the northern region of Ceará State is less susceptible to severe and persistent drought, such as the 2012-2018 event. Possible explanations are the intra-annual variability caused by proximity to the ocean. Another fact is that the higher latitudes are less impacted by ITCZ being positioned north of its climatological position. Therefore, in the central and southern regions, which contain 74% of the state's potential storage capacity, an extreme event can be more recurrent and water security can be compromised. Such analysis can be incorporated in drought plans to detect more exposed areas to drought. A limitation of the present approach is that it reflects meteorological drought exposure and does not consider water transfer between hydrographic regions, which may cause different levels of drought risk depending on where the supply is provided.
Classical approaches such as univariate analysis underestimates events frequency, while others cannot be associated with impacts if analyzed at inappropriate scales. This paper has argued that simultaneously considering drought duration and severity at a useful scale improves risk assessment. The return period, calculated as 240 years, reinforces the importance of proper documentation of the capabilities developed during this event once the same generation of decision-makers is not expected to face such an extreme event. The presented framework has shown that hydrographic region scale is adequate to couple drought impacts with the awareness given by bivariate analysis and can be replicated in drought plans for other regions. Copula functions were vital to jointly model drought characteristics, as other models cannot cope with their asymmetric behavior. Further investigation should analyze the scale that best represents specific impacts such as reservoir operation, water transfer between regions, and urban and agricultural supplies. Also, efforts should be made to understand the influence of events with different expected intervals with potential impacts on reservoirs levels, streamflow volumes and ecosystem thresholds. Several changes need to be considered in order to mitigate drought impacts; transforming statistical information into useful information for decision-makers is one of them.