Comparative evaluation of the dynamics of terrestrial water storage and drought incidences using multiple data sources: Tana sub-basin, Ethiopia

Ethiopia is referred to as the ‘water tower of Africa’ (MoWR 2011; Khadim et al. 2020; Woldegebriel et al. 2022) since many of its rivers originate from the central plateau and flow toward neighboring countries. The Blue Nile River is one of the major tributaries of the Nile River, which originates from Lake Tana (the largest lake in Ethiopia) and flows through Ethiopia and Sudan before joining the White Nile at Khartoum to form the Nile River.

The Tana sub-basin, located in the Upper Blue Nile Basin, has ample reserves of both groundwater and surface water resources (Ayenew et al. 2008; Abiy & Melesse 2017; Kindie et al. 2019). However, the sustainability of these resources has become questionable due to the increasing demand for agricultural, domestic, and industrial uses (Abiy et al. 2016a, 2016b; Abera et al. 2021). Overexploitation of these resources without considering sustainability may worsen the effects of human activities and climate extremes such as drought (Abiy et al. 2016a, 2016b; Nigate et al. 2016; Wakode et al. 2018; Abera et al. 2021). The stresses on groundwater and surface water resources are interrelated, and their depletion can negatively impact soil moisture as well (Abiy et al. 2016a, 2016b; Mamo et al. 2016; Walker et al. 2019; Tigabu et al. 2020). Therefore, proper management of terrestrial water storage (TWS), which integrates mainly surface water, groundwater, and soil moisture, is indispensable for sustainable development and proper utilization of water resources of the sub-basin. Improved water resources management requires monitoring systems that timely and accurately estimating the status of total or terrestrial water sources (Strassberg et al. 2009). Conventional methods for monitoring water resources and climate extremes (drought) involve taking hydrological and meteorological measurements at specific points, which may not be available or adequate in developing countries. Moreover, these methods are time consuming and expensive for large regions (Strassberg et al. 2009; Ochsner et al. 2013; Bi et al. 2014; Sun et al. 2021). Therefore, utilizing remotely sensed (RS) data and integrating it with traditional ground-based methods is essential to gaining a thorough understanding of the quantitative fluctuations of water resources in a given area.

Gravity Recovery and Climate Experiment (GRACE) (March 2002 to October 2017) and its successor GRACE-FO (GRACE-follow on) (since May 2018 to present) provide high-accuracy monthly gravitational field measurements of the water resources (Wang & Zhang 2017). Hereafter, GRACE is used to designate both GRACE and its follower GRACE-FO for convenience. GRACE has two identical satellites orbiting together at 220 km apart in a near-circular orbit and are connected by a highly accurate K-band microwave used to measure the distance variations between them due to gravity variations, including an accelerometer to remove nongravitational forces such as atmospheric effects (Tapley et al. 2004). GRACE is the recent remote sensing technology that paves a new way to assess the dynamics of terrestrial water storage anomalies (TWSAs). TWSA involves surface water storage anomaly (SWSA), snow water equivalent storage anomaly (SWESA), soil moisture storage anomaly (SMSA), and groundwater storage anomaly (GWSA) (Sun et al. 2021; Zhang et al. 2022).

Global Land Data Assimilation System (GLDAS) and land surface models (LSMs) are also used to estimate TWSA, which can be employed to compare with the GRACE-derived TWSA (Syed et al. 2008; Mohamed 2020; Liu et al. 2023). Numerous researchers have employed GRACE and GLDAS for different applications including assessing the dynamics of TWSA (Seyoum 2018; Xiong et al. 2021; Djessou et al. 2022), monitoring drought (Chen et al. 2018; Elameen et al. 2023; Getahun et al. 2023), and conducting groundwater change studies (Barbosa et al. 2022; Rusli et al. 2023) are a few examples.

Drought refers to a sustained TWS deficit caused mainly by climate variability, which severely influences socioeconomic development and environmental systems (Bayissa et al. 2018; Nigatu et al. 2021b; Getahun et al. 2023). Drought is becoming more frequent and intense in many parts of the world, including Ethiopia (Enyew et al. 2014; Nigatu et al. 2021b; Elameen et al. 2023; Getahun et al. 2023). Recurrent droughts occur in Ethiopia every 8–10 years or even 3–5 years, especially in the north (Enyew et al. 2014; Getahun et al. 2023). The Tana sub-basin in the Upper Blue Nile region has been the subject of several drought studies (Awange et al. 2014; Enyew et al. 2014; Legese et al. 2016; Bayissa et al. 2018; Seyoum 2018). Nevertheless, these studies have relied on traditional hydrological modeling that uses limited in situ meteorological and hydrological data. The use of RS data for investigating both TWS dynamics and drought monitoring is very limited in the region. Investigating the long-term TWS variations using multiple data sources and drought events using weighted water storage deficit index (WWSDI) is the first of its kind in Tana sub-basin. This study aims to (i) evaluate the long-term variations of TWSA employing different sources of TWSA, (ii) identify drought events and their interaction with TWSA dynamics, and (iii) test the performance of WWSDI in identifying drought events.

The GRACE dataset has four processing (using different algorithms) and archive centers, namely, Center for Space Research (CSR), Jet Propulsion Laboratory (JPL), Goddard Space Flight Center (GSFC), and GeoForschungs Zentrum Potsdam (GFZ) (Wahr et al. 1998; Landerer & Swenson 2012; Khorrami & Gunduz 2021). These GRACE satellite datasets have two forms, mass concentration (Mascon) and spherical harmonics (SH) solutions. The GRACE Mascon solutions have a lesser leakage error than the SH (Save et al. 2016; Scanlon et al. 2016; Wei et al. 2022). In addition, Mascon solutions do not need postprocessing, filtering, and have less dependency on scale factors (Save et al. 2016; Li et al. 2019). Thus, the GRACE Mascon solutions retrieved from CSR (http://www2.csr.utexas.edu/grace), GSFC (https://earth.gsfc.nasa.gov/geo/data/grace-mascons), and JPL (https://podaac-tools.jpl.nasa.gov/drive/files/GeodeticsGravity/tellus/L3/mascon) were employed.

GLDAS consists of LSMs using satellite and in situ observations to estimate disintegrated TWS components such as SWES, SWS (summation of canopy water storage and surface runoff (Qs)), and SMS (Sun et al. 2021; Xiong et al. 2021). GLDAS LSMs comprise the community land model, the Mosaic, the variable infiltration capacity (VIC), the NOAH, and the catchment land surface (CLSM) models (Li et al. 2019). In the present study, three GLDAS version 2.1 LSMs (Noah, VIC, and CLSM) were selected to calculate the average TWSA and compared with GLDAS 2.2-derived and GRACE-derived TWSA.

The latest version of GLDAS is GLDAS-2.2, which renders integrated TWS including groundwater storage (GWS) daily and 0.25-degree spatial resolution data from February 2003 to the present. It is the product of ground-based and model combined data, and data assimilation with CSR-GRACE/GRACE-FO (Li et al. 2019; Hasan et al. 2021). The CSR mascon solution was selected in GLDAS 2.2 data assimilation because of its finer spatial resolution, no postprocessing requirement, and less leakage (Li et al. 2019). The data were retrieved from https://disc.gsfc.nasa.gov/datasets?keywords=GLDAS_CLSM025_DA1_D_2.2&page=1.

Unlike the earlier GLDAS versions and GRACE, GLDAS-2.2 is the best for providing finer spatiotemporal resolution of TWS and GWS products. It also avoids the data gaps in TWSA derived from GRACE/GRACE-FO created between the GRACE end and GRACE-FO launch. In the present study, anomalies for the 20-year study period (from 2003 to 2022) were computed using the same base time (2004–2009) as GRACE-derived TWSA. The GLDAS 2 documented by Rui et al. (2020) was inferred for further information for all the aforementioned GLDAS versions and products. The Climate Data Operators, Quantum Geographic Information System (QGIS), RStudio, and python were applied for the GRACE and GLDAS NetCDF data processing and analysis.

The satellite altimetry missions including Topex/Poseidon, jason2, jason3, Jason, gfo, Envisat, Saral, and sentinel3A have provided an accurate and continuous surface height of lakes and inland (Crétaux et al. 2011; Ni et al. 2017). Several researchers have used satellite altimetry to estimate the lakes' water level changes (Calmant et al. 2008; Zhang et al. 2011; Ni et al. 2017). In the present study, the satellite altimetry data derived from the aforementioned missions was downloaded from http://www.legos.obs-mip.fr/soa/hydrologie/hydroweb/StationsVirtuelles/SV_Lakes/Tana.html and used. These data were employed to surrogate the data gaps in the ground-based Lake Tana water level gauge station (LTLGS) installed near the Bahir Dar meteorological station (Figure 1) for TWSA comparisons. The Tana water level obtained from combining the in situ and satellite altimetry was converted to monthly anomalies by removing the mean of the same base time (2004–2009) as GRACE-derived TWSA data.

Available data such as meteorological data (precipitation. temperature, wind speed, sunshine hours, and relative humidity) and daily water level of Lake Tana (at LTLGS) were acquired from the National Meteorological Agency (NMA) and the Ministry of Water and Energy (MoWE), respectively, for verification of GRACE- and GLDAS-simulated TWSA and drought indices.

Unlike GLDAS and altimetry data, GRACE data have 18 missing data from February 2003 to June 2017. There are also 11 consecutive months' data gaps (from July 2017 to May 2018) due to the gaps between GRACE ending and GRACE-FO launching periods. Besides, the August and September GRACE data were missed in 2018. First, the nonconsecutive missing data from February 2003 to June 2017 were filled using the linear interpolation method in the pandas-python data analysis library. Second, by using the previously filled data, the 11 consecutive months of missing GRACE data were filled by applying the seasonal auto-regressive integrated moving average (SARIMAX) forecasting model. Python package (statsmodels.tsa.statespace.sarimax) was employed for the imputation process (Barbosa et al. 2022).

TWSA provides a better opportunity to detect climate extremes than the traditional drought or wetness measuring indices such as standard precipitation index (SPI) and standard potential evapotranspiration index (SPEI) (Xiao et al. 2023). In recent years, researchers have made use of the opportunity of TWSA derived from GRACE to estimate droughts using the terrestrial water storage deficit index (WSDI) (Hasan et al. 2021; Nigatu et al. 2021b; Getahun et al. 2023). However, WSDI uses only the integrated TWSA and does not consider its components, which cannot provide better drought/wetness index results (Elameen et al. 2023). To better estimate the drought/wetness index, a WWSDI, which considers the contribution of each component of TWSA (Wang et al. 2020; Elameen et al. 2023), was adopted in the present study.

SPI and SPEI are the widely used traditional methods for drought monitoring (Seyoum 2018; Lotfirad et al. 2022; Getahun et al. 2023). SPI and SPEI are used to identify meteorological droughts for the short timescales, and hydrological and agricultural droughts for the medium and higher timescales (3 months and greater timescales) (Cammalleri et al. 2019; Getahun et al. 2023). These were employed here to compare the drought indices derived from TWSA. In the computation of SPI, monthly rainfall data were used as input, whereas in SPEI computation, monthly climatic-water balance, which is the rainfall after removal of potential evapotranspiration (PET), was employed. PET can be calculated using different available methods such as Hargreaves, Thronthwaite, Penman–Monteith, and FAO method procedures based on the available datasets (Beguería & Vicente-Serrano 2017). In the present study, Hargreaves and FAO methods were applied using SPI and SPEI functions in the SPEI RStudio package. The detailed procedures to calculate SPI and SPEI using spi and spei functions in the SPEI RStudio package are described in Beguería & Vicente-Serrano (2023).

The terrestrial water storage anomaly (TWSA) derived from CSR (CSR-TWSA), JPL (JPL-TWSA), and GSFC (GSFC-TWSA) with missed values are depicted in Figure S1. Missed data in each GRACE-TWSA source were imputed, and then the GRACE-TWSA, the ensemble mean of the three sources of GRACE, data were computed (Figure S2). GRACE-TWSA data are strongly correlated (r = 0.99) with CSR-TWSA, JPL-TWSA, and GSFC-TWSA data. Nevertheless, GRACE-TWSA is differently associated (NSE = 0.97, RMSE = 2.11), (NSE = 0.91, RMSE = 3.58 cm), and (NSE = 0.96, RMSE = 2.25 cm) with CSR-TWSA, JPL-TWSA, and GSFC-TWSA data, respectively. It can be perceived here that CSR-TWSA has the best correlation with the GRACE-TWSA in the Tana sub-basin. The previous study by Abiy & Melesse (2017) took only the TWSA derived from the CSR mascon solution and showed a plausible result in the Tana sub-basin. Nigatu et al. (2021b) also employed and confirmed that CSR-mascons have lesser uncertainty than the GRACE SH in the Nile River basin, which includes the present study area.

The average GLDAS-derived TWSA was estimated using GLDAS 2.1 products obtained from NOAH, VIC, and CLSM. The TWSA components, namely, SWS anomaly (SWSA), SMS anomaly (SMSA), and SWESA, were derived from the mean of these three GLDAS LSMs. Rivers, lakes, and reservoirs are included in the SWSA but are not specifically modeled in GLDAS models. On the other hand, Qs anomaly (QsA) can substitute for SWSA (Bhanja et al. 2016; Shamsudduha & Taylor 2020), or QsA and CWSA can be summarized (Sun et al. 2021; Xiong et al. 2021). Here, the total of QsA and CWSA (CWSA + QsA) was used. Less water equivalent depth (ranging from −0.014 to 0.038 cm) was seen in the result, which was not included in the GLDAS-derived TWSA computation method (Figure S3).

In general, the dominant components of TWSA derived from GLDAS 2.1 and GLDAS 2.2 are SMSA and GWSA, respectively. According to Wang et al. (2023), reducing the significance of SWSA in GRACE-TWSA decomposition by completely substituting it with the sum of QsA and CWSA might result in inaccurate findings. Twenty percent of the research area is made up of Lake Tana, and even in cases where reservoirs only make up one 1% of the watershed, their contribution to SWSA may not be disregarded (Longuevergne et al. 2013; Abiy & Melesse 2017). However, this study is not intended to address the process of translating Lake Tana level data to the real SWSA of the entire sub-basin (Djessou et al. 2022; Wang et al. 2023).

The GRACE-TWSA exhibits a good correlation (r = 0.88, NSE = 0.7, RMSE = 8.59 cm) with a month lag LTWSA and (r = 0.87, NSE = 0.73, RMSE = 7.06 cm) with 2-month lag PA. On the other hand, the CLSM_CSR-TWSA has the goodness of fit (r = 0.84, NSE = 0.57, RMSE = 10.15 cm) with a month lag LTWSA and (r = 0.88, NSE = 0.67, RMSE = 7.76 cm) with a 2-month lag PA. The GLDAS-TWSA has also a good relation with a month lag LTWSA (r = 0.64, NSE = 0.09, RMSE = 14.83 cm) and with 2-month lag PA (r = 0.79, NSE = 0.5, RMSE = 9.6 cm). The goodness-of-fit results revealed that GLDAS-TWSA has the least correlation with in situ LTWSA and PA. This is mainly attributed to the lack of the GWSA component. Therefore, it could not be taken for further analysis. The GRACE-TWSA has a bit better goodness of fit than CLSM_CSR-TWSA. This is because GRACE-TWSA has a lesser amplitude than that of CLSM_CSR-TWSA (Figure 6). Nevertheless, a lesser amplitude of GRACE-TWSA may be induced due to leakage errors in GRACE-TWSA (Wei et al. 2022). Moreover, unlike GRACE-TWSA, CLSM_CSR-TWSA has no missed values, and thus, no errors were induced during the missing data-filling process (Elameen et al. 2023). Therefore, CLSM_CSR-TWSA was selected for drought analysis.

The dominant components of the CLSM_CSR-TWSA are GWSA and SMSA as mentioned earlier. The advanced drought/wetness index (WWSDI) was calculated using these two components following the procedures suggested by Elameen et al. (2023). The CLSM-CSR-TWSA decomposition revealed a trend, seasonal/cyclic, and residual components, which contribute 56, 12, and 32%, respectively. These components were also used to calculate drought/wetness indices to compare with WWSDI (Khorrami & Gunduz 2021).

The WWSDI was compared with other commonly used traditional drought/wetness indices such as SPI, SPEI, and WSDI (Figure 7). In addition, the water storage deficit indices derived from the deseasoned (WSDIds), detrended (WSDIdt), and residual (WSDIres) values of TWSA were also employed to compare with the WWSDI as suggested by Khorrami & Gunduz (2021).

The WWSDI has a very high correlation with the GWSA-derived drought index (WSDIgws) (r = 0.997, NSE = 0.99, RMSE = 0.08 cm) (Table 2) since it is the largest component of the TWSA. Similarly, the drought magnitude (D_m) and intensity (D_int) of WWSDI throughout the study period (February 2003 to September 2022) have approximately the same as that of WSDIgws. For example, in the prolonged drought period from November 2004 to December 2006, D_m and D_int identified by WWSDI were 27.093 and 1.042, respectively (Table 1), whereas D_m and D_int identified by WSDIgws were 26.796 and 1.031, respectively. Second, WWSDI showed a strong correlation (r = 0.99, NSE = 0.98, RMSE = 0.13 cm) with WSDI, which was validated and widely used drought/wetness index around the globe (Sinha et al. 2017; Cammalleri et al. 2019; Nigatu et al. 2021b; Getahun et al. 2023).

The third-best correlation (r = 0.98, NSE = 0.95, RMSE = 0.21 cm) (Table 2) was observed with WSDIds. This was linked mainly to the least contribution (12%) of seasonal/cyclic component to TWSA, and thus, removing it was not sensitive for drought dynamics derived from WSDIds. WWSDI has also a fair correlation (r = 0.565, NSE = 0.13, RMSE = 0.93 cm) (Table 2) with the drought/wetness index derived from LTWSA (WSDIltwsa), a bit better than WSDI, which revealed the enhanced performance of WWSDI over WSDI (Table 2). In the present study area, however, WWSDI and WSDI are strongly correlated and can alternatively be employed to identify drought indices. This is mainly associated with GWSA accounting for a higher weight component ratio (CR = 0.73), whereas SMSA accounted for a lesser component ratio (CR = 0.27) of the TWSA.

It is interesting to note that the WWSDI) is considered to be a superior drought/wetness index compared to the WSDI. Unlike the WSDI, the WWSDI takes into account the contribution of TWSA components. However, in the present study, both WWSDI and WSDI showed similar results in identifying droughts, including onset, duration, intensity, and magnitude. This can be attributed to the fact that the GWSA component has a dominant contribution (73%) to TWSA.

WWSDI and WSDI have the best correlation with SPI9, SPI12, SPEI9, and SPEI12 within the sub-basin (Figure 8 and Table 3) and site-specific (Tables S1–S3) conditions. WWSDI and WSDI have relatively the highest association with SPI12 and/or SPEI12 of the Tana sub-basin and other site-specific except the Dangila meteorological station. The relatively highest correlation between WSDI and SPI12 was in agreement with the past study of drought using WSDI by Seyoum (2018). The rationale that WWSDI and WSDI have the highest correlation with longer timescales (9 and 12) was that GWSA is the dominant input of both WWSDI and WSDI, which has a good connection with slow hydrological responses (Cammalleri et al. 2019). More improved correlations have been observed between WWSDI/WSDI and SPEI indices than between WWSDI/WSDI and SPI indices (Table 3) due to PET effects. In general, the correlations between SPI indices and their corresponding SPEI indices are very high (ranging from r = 0.94 to r = 0.98). This portrays that the longer timescale hydrological droughts in the study area are mainly controlled by precipitation and less influenced by PET (Getahun et al. 2023).

In general, WWSDI and WSDI are reasonably related to SPI and SPEI indices, but have poor correlations (r ranged from 0.212 to 0. 435); details can be referred to Tables 1, S1, S2, and S3. The existence of these poor correlations was mainly attributed to the differences in hydrological responses between TWSA-derived indices (WWSDI and WSDI) and precipitation and PET-derived indices (SPI and SPEI) (Sinha et al. 2017).

In the present study, the long-term variations of GRACE-TWSA, CLSM-CSR-TWSA, and GLDAS-TWSA were employed and analyzed for the study period ranging from 2003 to 2022 in the Tana sub-basin. The GRACE-TWSA and CLSM_CSR-TWSA exhibited a very high correlation (r = 0.93) with each other. The TWSA obtained from the three aforementioned sources has also a strong temporal correlation with hydrological (Tana Lake water storage anomaly) and meteorological (precipitation anomaly) data. Unlike the GLDAS-TWSA, GRACE-TWSA and CLSM_CSR-TWSA showed a significant (p < 0.05) increasing (0.054 cm/month and 0.604 cm/year) and (0.042 cm/month and 0.431 cm/year) trend, respectively. CLSM_CSR-TWSA was chosen for drought analysis. The net annual CLSM-CSR-TWSA showed a terrestrial water storage deficit in 2005, 2006, and 2009 due to the climatic extremes (drought events) that occurred during those years. The drought events identified in those years and other years of a few months using WWSDI ranged from normal to extreme, and strongly influenced the dynamics of TWSA. However, the overall long-term trend analysis showed a significant (p < 0.05) increasing trend and 50.68 cm net equivalent water height gain in the study period. The WWSDI was compared with drought indices such as WSDIgws, WSDIds, WSDIltwsa, and traditionally used indices (including WSDI, SPI, and SPEI). In this particular study region, WWSDI has a very high correlation (r = 0.99) with the commonly used drought/wetness index (WSDI) and poor correlation with SPI and SPEI at different timescales. The WWSDI is an advanced drought/wetness index and is better than other indices since it considers all the TWS components. Although the overall TWSA showed a significant increasing trend, this study recommends different drought adaptation techniques including afforestation and water conservation to lessen water losses concomitant to droughts. This study may also be used as the basis for future researches focusing on spatiotemporal water storage and drought dynamics in the studied region. This study, in general, suggests that employing GRACE- and GLDAS-derived TWSA is vital for analyzing the long-term dynamics of total water storage and drought monitoring for sustainable water resources management, especially in data-scarce regions.

We are grateful to Arba Minch University, Water Resources Research Center who partially financed this study under a small grants scheme with a project code of GOV/AMU/TH28/AWTI/WRRC/03/2014. We also thank the Ethiopian Ministry of Water and Energy for providing the hydrological data and the Ethiopian National Meteorological Agency (ENMA) for providing the weather data. Finally, we greatly acknowledge Amhara Regional State Water Bureau and design offices for providing the field and historical water level data.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.