ASSESSMENT OF HYDROLOGICAL MASS LOSSES IN THE NORTHEAST ATLANTIC EASTERN HYDROGRAPHIC REGION, BRAZIL

Freshwater monitoring globally is fundamental to support decision-making. However, long-term hydrological data for some regions are lacking due to limited of observational networks. Remote sensing products come to aggregate the in situ observations and overcome obstacles about data availability. This study assesses the hydrological mass losses in the Northeast Atlantic Eastern Hydrographic Region (NAEHR) in Brazil through temporal data sources. GRACE RL06 Mascon solutions, reservoirs volume, vegetation index and rainfall stations are used. The results confirm the cohesion between the TWS variations with water volume and NDVI, obtaining a strong correlation of 0.78 and 0.77 respectively. The Standardized Precipitation Index scales (12, 24 and 36 months) showed a moderate and strong correlation with the TWS of 0.57, 0.73 and 0.72 respectively and identified the last regional severe drought from 2012 to 2017. The NEAHR is located almost entirely in the Brazilian semiarid region, comprises about 24 million inhabitants, it is high vulnerable to drought, thus justifying the importance of monitoring its water resources availability.


Introduction
It The hydrological drought events have been frequently occurring during the last decade. It has been possible to observe water demand influences among human consumption, agricultural irrigation and ecological systems stability (Bai and Liu, 2018;Jassey et al., 2018;Gentirana, 2016;Van Loon, 2015). Although, drought could be often attribute to a natural phenomenon, some studies indicated human-induced impacts (Cook et al., 2009;Van Dijk et al., 2013). In recent decades, droughts increased frequency and intensity all over the planet and climate changes could be correlated to it (Ndehedehe and Ferreira, 2020a;Zhou et al., 2011). The precipitation scarcity contributes to the drought phenomenon, consequently the low evapotranspiration and low soil moisture (Mpelasoka et al., 2018;Feng et al., 2019;Yao et al., 2020). In Northeast of Brazil, this phenomenon occurs mainly in the semi-arid region due to uneven precipitation in space and time (Cunha et al., 2015;Brito et al., 2017;Ndehedehe and Ferreira, 2020b). Brazil has an extensive water reserve, on average 12% of the total existing in the world, surpassing even some continents, such as Oceania, which holds 6%, Europe with 7%, and Africa with 9% (UNESCO, 2009). However, in some Brazilian Hydrographic Regions (HR), drinking water reserves are extremely low, such as the Northeast Atlantic Eastern HR (NAEHR), selected as study case and which has 0.2% reserve, the lowest water availability compared to other existing HR in the country (ANA, 2013). The NAEHR is located in the extreme northeast of Brazil, which runs from Piauí to Alagoas states, has around 286x10³ km² and 87% of this area is located in the Brazilian semiarid region.
The drought is normally classified into four categories, namely as meteorological, hydrological, agricultural and socioeconomic (Wilhite and Glantz, 1985). The main drought causes are due to: a) rainfall deficiency in particularly the distribution and intensity, b) storage and c) demand and the existing water use (Zarch et al., 2015). Therefore, to mitigate climate change impacts, it is necessary to have a better understanding about the drought intensity (duration and spatial extent), the seasonal atmospheric variability and the drought severity trends, among other issues (Shabbar and Skinner, 2004).
There are several indexes for the drought identification and classification, among them, the Standardized Precipitation Index (SPI) (McKee et al., 1993). This index makes use of statistical temporal rainfall data analysis and has been widely applied to rainfall assessment, monitoring and prediction (Spinoni et al., 2014) e.g., in regions such as Africa (Cheo et al., 2013;Dutra et al., 2013), China (Feng et al., 2013), India (Dhakar et al., 2013) and United States (Ford and Labosier, 2014). This versatility allows SPI to provide insights for the agricultural production, water resources and groundwater management (Mishra and Singh, 2010).
Another parameter related to drought detection is the availability of the Total Water Storage (TWS) which is the virtual sum of all water stored above and below the Earth's surface including water in lakes, rivers, human-made reservoirs, wetlands, soil, and groundwater reservoirs (Cazenave and Chen, 2010). Awange et al. (2014) applied the inter-annual and seasonal TWS variability over Ethiopia. Ahmed et al. (2014) used the Gravity Recovery and Climate Change (GRACE) data to monitor natural and man-made induced variations in water availability across Africa. Ndehedehe et al. (2016) presented a study case in West Africa during the period 2002 to 2014 using the GRACE's monthly TWS. In Australia, Chen et al. (2016) studied long-term groundwater storage changes in Victoria, using GRACE and in situ observations. Molodtsova et al. (2016) assessed the GRACE potential for flooding mitigation in the United States of America. In China, Jiang et al. (2017) studied seasonal drought combining crustal deformations from GPS data and TWS. In Asia, Xiang et al. (2016) studied the TWS changes for the Tibetan plateau and adjacent areas. Sinha et al. (2017) used the GRACE observations to characterize the drought in India between April 2002 and January 2015. In Brazil, Sun et al. (2016) used GRACE data, from 2002 to 2015, to assess the ongoing drought in the São Francisco River watershed. Montecino et al. (2016) demonstrated another GRACE applicability for South America with groundwater focus considering the northern region of Chile. Rosenhaim et al. (2018) studied the behavior of TWS in the Northeast Atlantic Eastern Hydrographic Region using GRACE, TRMM and in situ observations for the period 2002 to 2015, thus identifying a loss of approximately 23.710hm³ between the extreme years studied.
The use of remote sensing can complement hydrographic studies such as the Moderate Resolution Imaging Spectroradiometer (MODIS) on board of the Aqua and Terra satellites. In this case, the drought quantification, take advantage of spectral indices and water balance simulations (Cunha et al., 2015). The MODIS products used for terrestrial applications, aimed to monitoring the planet's vegetation cover, able to identifying changes, either resulting from climate or weather changes (Justice et al., 2002). Deus and Gloaguen (2013) combined MODIS data and the Modified Normalized Difference Water Index (MNDWI) to identify the water changes over Lake Manyara, Tanzania. Son et al. (2012) used MODIS to monitor agricultural drought in the Lower Mekong Basin in Southeast Asia that covers Thailand, Laos, Cambodia and Vietnam. Gu et al. (2007) used the MODIS sensor to assess drought by physical indexes in the central United States of America. Likewise, Andam-Akorful et al. (2017) used the vegetation index retrieved from the Advanced Very High Resolution Radiometer (AVHRR) to assess the human influence on the water resources of West Africa. In Brazil, Anderson et al. (2010) used the MODIS sensor products to study the Amazon forest crowns drought, whereas Morton and DeFries (2005) used MODIS to assess the annual Amazon deforestation.
The hydrological changes reported previously provide methodological support for patterns detection and trends analyses around the world. Thus, a new contribution for NAEHR, as study case combining in situ and remotely sensed products are used to understand the water mass variability, which the objectives to: (i) identify the spatiotemporal pattern of the total water storage changes and its possible trend for the next years using GRACE Mascons solutions; (ii) analyze the temporal reservoirs volume; (iii) extract and analyze the Normalized Difference Vegetation Index (NDVI) from MODIS; (iv) detect the extreme drought events using SPI at different time scales and; (v) perform the statistical correlation between: a) TWS and reservoir volume, and b) TWS and NDVI.

Northeast Atlantic Eastern Hydrographic Region (NAEHR)
The (NAEHR) (see Fig. 1b) is one of the 12 hydrographic regions in Brazil (see Fig. 1a) that occupies 3.4% of the national territory with an area of approximately 286,800 km², covering six states: Piauí, Ceará, Rio Grande do Norte, Paraíba, Pernambuco and Alagoas (ANA, 2015). Its location (almost all) belongs to the Brazilian Semi-Arid Region, which the lowest water availability in Brazil, with prolonged droughts periods and high temperatures throughout the year. In addition, it has a great economic and social importance for the region, covering about 740 municipalities and five state capitals (Ceará, Rio Grande do Norte, Paraíba, Pernambuco and Alagoas). The region has its greatest water demand considering urban supply and irrigation (ANA, 2018).
In addition, it comprises 24 million inhabitants, thus representing about 12% of the country's population (ANA, 2015).
The study area includes fragments of the Atlantic forest, Caatinga, Cerrado, and coastal biomes. The Caatinga have been devastated over the decades by livestock that invaded the hinterlands and the forest zone deforested by the sugarcane culture implantation. Even today, plant extraction, mainly for the exploitation of timber potential, represents one of the activities affecting the environment (ANA, 2015).
According to the Water Resources Conjuncture in Brazil (ANA, 2015), the annual NAEHR rainfall average is 1.052 mm, below the national average, 1.761 mm. The water surface availability, considering the flow regularized by the region's reservoirs is 91.5 m 3 /s, which corresponds to 0.1% of the country's surface availability (91.071 m 3 /s). The average flow of the Hydrographic Region (HR) is 774 m 3 /s, corresponding to 0.43% of the national average flow (179.516 m 3 /s).
The HR has low groundwater storage capacity, due to the low aquifers recharge and rainfall (Hu et al., 2017). In this region the crystalline rocks consist more than 90%. Thus, there is a certain groundwater predominance with high salt content. Normally the sedimentary regions, with better water quality are in the coastal areas (Cirilo, 2008).

Dataset
All data (multispectral images from the MODIS sensor, GRACE RL06 Mascon solutions, total reservoir volume and precipitation) applied to study the hydrological mass changes at NAEHR were obtained free of charge. The sections 2.2.1 to 2.2.4 showed all dataset baselines

Reservoir Volume
The reservoirs volumes are obtained by the Reservoir Monitoring System of the Brazilian National Water Agency (ANA) (https://www.ana.gov.br/sar0/MedicaoSin). There are 346 reservoirs existing in the study area distributed among six states (Alagoas, Ceará, Rio Grande do Norte, Paraíba, Pernambuco and Piauí). Table 2 shows the states distribution, capacity and participation of each state. The data selected consists the daily volumes of the total reservoirs from April 2002 to June 2017.

Moderate Resolution Imaging Spectroradiometer (MODIS)
The MODIS sensor board on the TERRA satellite monitors changes on the land surface among others applications (https://search.earthdata.nasa.gov). The product MOD13 selected provides consistent information on vegetation conditions, with a temporal resolution of 16 days and spatial resolution of 1000 m. To cover the NAEHR are necessary three scenes at the end of the mosaicking process. It were select 78 scenes described in section 2.3.3.

Rainfall
The monthly-accumulated precipitation data refer to a 30 years study period from August 1989 to July 2019 according to the 26 rainfall stations provided by the Brazilian National Meteorological Institute (INMET) -shown in Figure 1b.

Method
This section discusses the methods and processes used in the research divided into four stages: a) data; b) pre-processing; c) processing and d) analysis, as shown in Figure 2.

Spatio-temporal patterns
In order to find the dominant spatio-temporal patterns, the Empirical Orthogonal Functions (EOF) method or also known as PCA (Principal Component Analysis) was applied. The EOFs have been used successfully in the analysis and modeling of temporal variations in mass distribution obtained by the GRACE satellites (Crossley et al., 2004;Li et al., 2016;Peng et al., 2017, among others). This method of statistical decomposition seeks to reduce the size of a data set by finding some linear orthogonal combinations of the original variables with the greatest variance and preserving the main variation of the original variables. In summary, EOF analysis allows to represent the TWS time series as: where TWS is the original time series as a function of time (t) and space (ϕ,λ) of the main factors responsible for the temporal variations of TWS while the Principal Components (PC) shows how the amplitudes of each EOF vary with time.
The GRACE Mascons time series (monthly and annual averages), from April 2002 to June 2017, are also used to extract the linear trend by the ordinary least squares (OLS) method (Gemael, 1994). Thus, it was possible to reach the objective (i) described in the introduction.

Reservoirs volume analysis
The Brazilian national water agency defines the reservoirs operating rules in Brazil and present monthly monitoring bulletins. The daily values sum for the 346 reservoirs, obtained the useful monthly volume time series, also the annual averages and monthly maximum and minimum. The dispersion graphs were prepared together with the linear regression in Excel software to obtain its analytical equation and trend line, regarding the objective (ii).

Normalized Difference Vegetation Index -NDVI
The MODIS images have high temporal resolution, being georeferenced and corrected with atmospheric effects, contributing to a better understanding of global systems Justice et al., 2002). The images in the HDF format, in the pre-processing stage are converted to GEOtif format using the MODIS Reprojection Tool software. The date choice is according to the months matching with the highest and lowest value of TWS from mascons GRACE, for each year of the time series, resulting in the months selected on table 3. With the pre-processed images inserted in the software QGIS, it was possible to select the study area and process the NDVI. This index used to represent the vegetation cover in response to the region's water variation. The NDVI range between -1 and +1 using bands of red and near infrared (NIR) based on the extent of reflectance surface (Carlson and Ripley, 1997). With the calculated index, the average of the NDVI values for each month is used to attend the objective (iii).

Standardized Precipitation Index -SPI
It was selected 26 rainfall in-situ stations (Fig. 1b)  With the re-structured dataset, it was possible to apply the standardized precipitation index -SPI, which consists of adjusting the gamma probability density function to a given frequency distribution of the total precipitation for a season (more details in Zarch et al., 2015). The SPI was evaluated according to 3-time scales, 12, 24 and 36 months with the idea to verify the behavior of shorter droughts, such as meteorological ones that are identified in smaller SPI time scales, to more lasting droughts such as hydrological ones in larger SPI time scales, important guidelines for water resource management (Bonaccorso et al., 2003). Awange et al. (2016) reported that short-timescale droughts are related to agricultural drought and river discharge in headwater areas, medium-timescale droughts reflect the reservoir storages levels and rivers, and relate to groundwater storage are the long-timescale droughts. There are no general recommendations on which scale is appropriate, it is always depending on the user experience and applications involved (Cheval, 2015).
The SPI Generator (https://drought.unl.edu/droughtmonitoring/SPI.aspx) software developed by the United States National Drought Mitigation Center (NDMC) was used to process the temporal series. The positive SPI values indicate precipitation above the median and negative values indicate precipitation below the median according to table 4, thus meeting objective (iv).

Spearman rank correlation
The Coefficient of Variation (CV) is calculated to identify homogeneity and heterogeneity of the variables (Snedecor and Cochran, 1980) (annual averages of the mascons GRACE, total volume of the reservoirs and average pixels of the NDVI), where values greater than 30% of CV were identified indicating that the variables are not parametric, so it is recommended to perform the Spearman rank correlation (Siegel, 1975). According to Schober et al. (2018), the Spearman rank correlation (ρ) does not require normally distributed data and can be used, analogous to Pearson coefficient varying from -1 to +1. The sign indicates the direction of the relationship between the variables, the negative coefficient indicates that the correlation is inverse, as one of the analyzed parameters increases the other decreases. The correlation is tested at a significance level of 1% bilateral analysis, thus indicating whether the variables have a statistically significant difference. It is adopted, as a null hypothesis assuming that the statistical result of the correlation was obtained by coincidence due to the probabilistic fluctuations of the events, and their alternative hypothesis considers a real similarity of the correlation. This stage allows extracting information to meet the objective (v).

Results and Discussions
All correlation analyzes discussed in this section obtained a p-value = 0.00, that is, less than 1%, thus accepting the alternative hypothesis that states that there is a real similarity in the correlation between the variables studied. Figure 3 shows the annual averages of the TWS in red and the annual average value of the useful volume of the reservoirs in blue, as well as the linear regression (trend lines) and equations for both samples. According to the results obtained for the regression equations, the slope coefficients are negative, indicating a successive decline from 2011 to 2017, despite being different sample sets, both have similar behavior over time. The TWS and the reservoirs volume follows the same trend. The data were subjected to Spearman (σ) where it was possible to evaluate the behavior between these two non-parametric and independent variables, whose result was (ρ) = 0.78. According to figure 3 it is possible to verify large peaks detected for the years 2004, 2009 and 2011 that may be correlated to various atmospheric phenomena. For example, 2004 was considered a very rainy year in the first months, according to Alves et al. (2006) in January 2004, the northeast trade winds became more intense than the southeast trade winds, driving the Intertropical Convergence Zone (ITCZ). In 2009 it was considered the rainiest year of the decade, with the highest rain concentration in April and May, which corroborates with Alves et al. (2009) the ITCZ convection was reinforced by the favorable pulse of the Madden-Julian Oscillation (MJO) over the northern of the Northeast Region between late April and early May. In 2011, the La Niña phenomenon event had a configuration called Modoki, named by the second variability mode of the Sea Surface Temperature (SST) Tropical Pacific more details in Martins and Vasconcelos Junior (2017). The year 2010 was classified as a dry year, only in 2011 did the northeastern Brazilian region showed rainfalls above normal (Marengo et al., 2018). This event is associated with another major drought in the Northeast, highlighting the year 2012 as the beginning. The linear regression indicates a decline tendency over the years, resulting in the lack of water for the region as well as the volume loss in the reservoirs. Considering the same study case Rosenhaim et al. (2018) identified the variation in NAEHR's water availability in the first and second semesters of the years 2002 and 2015, where the variation in the first semester was 23.71km³ and for the second half was 19.25 km 3 , representing a considerable loss in the estimated volume of water for this region.  Figure 4 shows the comparison between TWS and NDVI following the linear regression of the samples that resulted in the linear trend. The NDVI average vary between 0.34 to 0.74, in one hand, values closer to 1.0 indicate a large amount of vegetation and active photosynthesis. On the other hand, the values closer to 0.0 indicate little or no chlorophyll activity. The variation in vegetation activity has been associated with climate change (Sun et al., 2015), where vegetation connects with soil, atmosphere and water, serving as one of the climate change indicators (Cramer et al., 2001). According to 26 rainfall stations data observed in this research, the annual average rainfall is 2297.70 mm, where the first annual semester has precipitation values above the average and the second semester with precipitation values below the annual, thus characterizing two distinct seasons over the years. The NDVI response corroborates whit the region's rainfall, however within the time series there are anomalous rainy and drought months that may interfere with the correlation between NDVI and TWS. The NDVI values corroborate with the rainfall regime and the TWS variations in the region, the highest NDVI values occurred between March and June, where the rainiest months of each year are considered. Figure 4 in one hand shows 2002, 2003, 2006 to 2011 the results for NDVI ranging between 0.69 and 0.74 and positive variations of TWS between +0.99 and + 17.4cm. On the other hand, months like January and February that are part of the first half of the year, presented low values for NDVI, as in 2003NDVI, as in , 2004NDVI, as in , 2006NDVI, as in , 2008NDVI, as in and 2009 to the fact that NDVI does not have an immediate response to like precipitation. The periods of low rainfall as November and December reflected low NDVI values for the months of January and February. It is possible to observe that from 2012 until 2017 the negative values of the TWS are recurrent, where the NDVI average values were between 0.4 and 0.6. Thus being able to perceive a progressive NDVI decline from low peak values between 2014 and 2017. With this decay it is highlighted a waring alert of environmental impacts such as surface erosion (Singh et al., 2004;Asis and Omasa, 2007;Karaburu, 2010) and desertification (Piao et al., 2006;Feng et al., 2015;Tomasella et al., 2018). One of the reasons for the decline in these values (Figure 4) is the lack of rainfall, since water is a fundamental element for the vegetation health, whose terrestrial ecosystem hydrological control is exercised by the soil water balance (Rodriguez-Iturbe et al., 2007). According to Omute et al. (2012), the NDVI is able to describe the variability of water levels in response to weather patterns. The combination of TWS, precipitation and NDVI are very sensitive to vegetation drought in semiarid areas (Zhao et al., 2017). Considering the Spearman (ρ) was possible to evaluate the behavior between these two nonparametric and independent variables, whose result was (ρ) = 0.77. Corroborating with the results (i) Geruo et al. (2015), indicating that NDVI is influenced by TWS, and by water control associated with vegetation growth and (ii) Ma et al. (2012) emphasizing the vegetation growth reduction related to drought-induced water stress.  Figure 5 shows the three temporal variations of the SPI: figure 5 (a) refers to the 12 months SPI, which is directly associated with the water lack in flow water reserves and groundwater levels, which correspond to long time scales. Awangue et al. (2016) showed a 12-month drought probability ranging from 19% to 24% in the region under study, with a 6% to 10% extreme droughts probability of occurrence. The calculated SPI-12 identified a weak drought for the year 2010 corroborating with Marengo et al. (2018) that defined the year 2010 as a dry year for the northeast region. The beginning of the 2012 drought was associated with the La Niña phenomenon (Rodrigues and McPhaden, 2014), according to Marengo et al. (2018), the La Niña phenomenon only interfered in 2015-2016 aggravating the drought situation, where it is possible to identify such behaviors in figure 5 (a). The TWS and SPI-12 had the same behavior throughout 2002-2017 however, diverging in some points as the first half of 2007, 2010 and 2012. This is probably to the fact that SPI used at least 30 years precipitation time series and the TWS monthly variation of the available data. At the first semester of 2007, 2010 and 2012 low levels of precipitation were recorded whose SPI-12 classified them as dry months, however the positive precipitation variations occurred in April and May thus caused a positive TWS variation in June of those years.

Drought characterization using Standardized Precipitation Index
The Spearman correlation (ρ) helped to evaluate the behavior between the parametric and independent variables SPI-12 and TWS, whose result was moderate with (ρ) = 0.57. Figure 5 (b) refers to the SPI-24, being used to obtain the spatial distribution of drought and hydrological humidity (Bordi et al., 2009). Awange et al., (2016) showed the probability of drought for 24-month ranging from 16% to 22% in the study region, with a 7% to 12% probability of extreme droughts occurrence. The calculated SPI-24 was unable to identify the 2010 dry year, which according to Marengo et al. (2018)  Figure 5 (c) refers to SPI-36, since it is also a long-term scale, SPI-36 is associated with hydrological drought studies, which according to Abatan et al. (2017) are associated with persistent droughts. The index values on this 36-month time scale are higher, and at no point in the time series did the index identify severe drought. The Spearman correlation (ρ) = 0.72 represented a strong correlation between SPI-36 and TWS.
The three scales studied in this research are long-term, applied for hydrological droughts studies. The SPI makes the precipitation comparison of the consecutive months for each scale using the time series. For all cases the SPI-12-24-36 showed a moderate and strong correlation witch TWS. The scales of 24 and 36 months do not describe the drought over 2002-2019 in detail. The 6-year long drought that occurred between the years 2012 to 2017, was identified in all SPI scales used in the methodology, this was due to the duration of the drought event that exceeded 36 months, such scales are represented by Figure 5

TWS Space-Time Analysis
The EOF analysis identified that the first three modes are capable of recovering 98.8% of the TWS variance in the study region, more specifically PC1, PC2 and PC3 represent 82.2, 9.7, and 6.9% respectively. The rest of the PCs (<1.2%) were considered as noise, and then disregarded in this analysis. (82,2%), PC2 (9,7%) and PC3(6,9%) of TWS.
The PC1 represents the temporal variation of the annual, inter and intra-annual signals. The Lomb-Scargle periodogram of PC1 revealed that the predominant frequencies are approximately 1 cycle per year (cpy), 0.4 cpy and 2.2 cpy. EOF1 (see. Fig. 7a) presents the spatial distribution of the temporal patterns found in PC1. Furthermore, the annual cycle and the other cycles observed in PC1 showed greater variability in the western region, this variation increases progressively from the coast to the interior of the South American continent. These variations can be explained due to the significant annual phase (see. Fig. 8b) differences (up to 132 days) found in some parts of the study region, as well as the other intra and inter annual frequency signals. The EOF1 patterns resulted close to those found in the annual amplitude map (see. Fig. 8a) considering that places with the largest annual amplitude showed the highest EOF variations (in magnitude) and the ones with the least amplitude varied less.
The amplitude and phase of annual cycle were estimated based on the following functional model: where a 0 is the constant term, a 1 is the linear velocity, t is the time, 1 = 2 1 , 1 = 1 year, and S k and C k are the coefficients for computing the amplitude of seasonal signals 1 = √ 1 2 + 1 2 . The parameters of the functional model were estimated by least squares adjustment. EOF2 (see Fig. 7b) is associated with PC2 (see Fig. 7e), which represents the annual variability 1 cpy and semiannual 2 cpy of the TWS, this indicated in the periodogram. The spatial distribution of this variability indicates a gain in the north-west region, and a decrease in the south-east region of the study case. The EOF3 represents the spatial variation of the temporal variations associated with the annual cycle and the intra-annual variation of 2.27 cpy observed in the PC3 periodogram (see Figs. 7e and 7f). The largest variations are concentrated in the central/ southern-western region, and the smallest are concentrated in the northern-western of the study.

Conclusion
This study was carried out based on the MODIS sensor images (NDVI), rainfall data, useful reservoir volumes and Mascons GRACE observations, all combined to assess the hydrological mass behavior in the NAEHR, Brazil. Through the correlations and seasonal peaks identification of water availability, it was possible to demonstrate an alternative for the qualitative and quantitative assessment of water resources. The main results found indicate: i.
The Spearman's linear correlation (ρ) between the independent variables TWS with annual average value of the reservoirs (Figure 3) and the NDVI average (Figure 4), were 0.78 and 0.77 respectively, showing a strong correlation between them. This highlighted a direct correlation among TWS, reservoirs volume and NDVI, presenting outstanding parameters for water management in this region. ii.
The year 2009 obtained the largest reservoirs volume of water, with 16816.28hm³ and TWS = 10.24cm. The month of May presented the highest contribution for this variability, showing a volume of 18,331.59hm³ and TWS = 23.01 cm. The EOF's analysis enabled a better understanding of the spatial distribution and temporal changes of TWS in the study region and the respective frequencies associated with these variations. iii.
In general, the NDVI corroborates with the variations presented by TWS. Both variables showed the same behaviour throughout 2002-2017. The positive variations in TWS reflected on NDVI values close to 1.0, as occurred specifically in the first half of each year, thus indicating a dense vegetation cover that normally needs a greater hydric demand. On the other hand, NDVI also response to the dry period identified in the study, e.g., the results between 2012-2017 dropped in maximum values varying between 0.4 and 0.6.

iv.
The values obtained for each SPI time scale were able to identify the last drought period between the years 2012 and 2017, this index responded cohesively to the Mascons linear regression, which indicated the declining tendency considering the following years. The Obtaining Spearman correlation (ρ) showed of 0.57, 0.73 and 0.72 for the 12-month, 24-month and 36-month sales, respectively. Although the SPI-36 had the strongest correlation, the SPI-12 detailed better the drought temporal distribution.
v. The GRACE mission satellites come to assist also regional management of water resources in Brazil and its potential and dissemination for new study cases. The main advantage is the applicability for the inaccessible places where is not possible to collect in situ information. As disadvantages, the time series is still considered small for hydrological studies, since its operation started in 2002, which can be mitigated with the combination of other products as proposed in this work. Another limitation is the spatial resolution, which is not adequate when studying a small hydrographic basin or even a specific reservoir. The GRACE mission has its successor, Gravity Recovery and Climate Experiment Follow-On (GRACE-FO), operating since May 22, 2018. duration and severity of drought in the Semiarid Northeast Brazil region. International Journal of Climatology,v. 38,