Integrating GRACE products and land surface models to estimate changes in key components of terrestrial water storage in the Nile River Basin


 The Nile River Basin (NRB) is facing extreme pressure on its water resources due to an alarmingly increasing population that is extremely vulnerable in aspects of irrigation and hydropower. The NRB ascends itself to remotely sensed approaches with high resolution of spatial and temporal coverage as disparate to ground-based in-situ observations due to its size and limited access from basin countries. The Gravity Recovery and Climate Experiment (GRACE) allow a unique opportunity to investigate the changes in key components of Terrestrial Water Storage (TWS). Differences in tuning parameters and processing strategies result in GRACE TWS solutions with regionally specific variations and error patterns. We explored the spatiotemporal changes of the TWS time series, trend, uncertainties, and signal-to-noise ratio (SNR) among different GRACE TWS. We had also investigated the key terrestrial water storage components (surface water, soil moisture, and groundwater storage changes). The results show that the uncertainty of GRACE spherical harmonic (SH) solutions are higher than the mass concentration (mascon) over the NRB, and the Center for Space Research-mascons (CSR-M) noted the first best performance. Substantially, significant long-term (2003–2017) negative groundwater and soil moisture trend demonstrates a potential depletion over NRB. Despite an increase in precipitation and TWS time series, the rate of decline noted to increase rapidly from 2008, thus indicating the possibility of human-induced change ( e.g., for irrigation purposes). Thus, the result of this study provides a guiding principle for future studies in TWS change-related hydro-climatic change over NRB and similar basins.


1.
Introduction The Nile is the world's longest river. It has a drainage area of about 3.2 million km2 that is nearly 10% of the landmass of the African continent. It flows across 11 countries from South to North with approximately 35 degrees of latitude, which crosses highly diverse landscapes and climatic zones. The two main tributaries of the Nile River basin (NRB) are the White Nile and the Blue Nile. NRB is facing extreme pressure on its water resources due to alarmingly increasing population (Siam et al., 2017) that is extremely vulnerable in aspects of irrigation, agriculture (Abtew et al., 2014;Multsch et al., 2017), water supplies, and hydropower (Ahmed et al., 2019;Hasan et al. , 2018). Owing to its size, the Nile inclines itself to remotely sensed approaches with the high spatial and temporal coverage as the NRB disparate to ground-based in-situ observations. Therefore, understanding its surface, soil moisture, and groundwater storage are crucial to providing useful information for proper water resource utilization and climate change impact studies.
The Gravity Recovery and Climate Experiment (GRACE) offer a unique opportunity to quantify and investigate the Terrestrial Water Storage (TWS) variations at regional and global scales.
GRACE-derived TWS delivers vertically integrated water storage changes in all water-bearing layers. The integration, as indicated in Eq.1 (Ramillien et al., 2014), comprises four components.
These are surface water storage in rivers, lakes, and wetlands (∆ ), soil moisture storage (∆ ), ice and snow water storage (∆ ), and groundwater storage (∆ ). Thus, GRACE measurements have become a crucial hydrological tool for quantifying basin-scale TWS (Hu & Jiao, 2015). We used only three components (Castellazzi et al., 2016;Long et al., 2015) to assess spatiotemporal changes in TWS, assuming the contribution of ice/snow to be negligible over the NRB (Eq.1).
TWS represents the sum of all water stored above and below the Earth's surface. It is a significant component of the global hydrological cycle, counting water in lakes, rivers, wetlands, soil, humanmade, and groundwater reservoirs. TWS change controls the weather through a series of processes and feedback mechanisms that reflects the natural variability of climate (Yang et al., 2015). Thus, quantification and monitoring of TWS change are critical for characterizing changes in water resources and for improved prediction of regional-global water cycles and interactions with the Earth's climate system (Famiglietti, 2004).
GRACE's products are mainly in the form of traditional spherical harmonic (SH), and mass concentration (MASCONS) produced from different centers. These centers processed with different approaches, such as the variational equation integration approach, short-arc integral approach, energy conservation approach, celestial mechanics approach, and acceleration approach.
The major drawback of the GRACE SH solutions is the coarse spatial resolution. In another way, spatial filtering required to suppress the spatial noise "over degree and order 30" for GRACE SH coefficients that are dominated by the longitudinal stripping noise (Swenson & Wahr, 2006).
These drawbacks further degrade the spatial resolution of GRACE's estimated mass changes .
Thus, GRACE mascon solutions had released to overcome the drawbacks of the SH. In GRACE mascon solutions, mass change products are not affected by the longitudinal stripes and do not require smoothing or filtering, and gives improved spatial resolution Watkins et al., 2015). According to Wiese et al. (2016), the mascon solutions constrained by a priori information derived from near-global geophysical models to prevent striping and filtered by Coastline Resolution Improvement (CRI) to reduce spatial leakage from land to oceans.
Additionally, TWS of different GRACE solutions has studied at the global and basin level.
Ensemble prediction and inter-comparisons of the SH solutions were recommended by (Sakumura et al ., 2014) due to reduced noise in the gravity field solutions within the available models. At the global level (Scanlon et al., 2016);  studied the global trends and annual amplitudes by comparing the mascon with SH solutions and recommended mascon solutions for accurate surface-based gridded information used without further postprocessing. At basin level, (Shamsudduha et al., 2017) used JPL-M and GRACE Tellus products that highlighted the substantial uncertainty in the amplitude of Terrestrial Water Storage Anomaly (TWSA). None of those mentioned above studies considered a recently released GSFC mascon solution. Some river basins experienced discrepancies between GRACE solutions (Scanlon et al., 2018) that leads to significant inconsistencies in TWSA trend estimations. Hence, evaluating those discrepancies is urgent for understanding uncertainties in monitoring water resources variations based on basin scales.

GRACE and GLDAS (Global Land Data Assimilation System) have provided useful observations
for studying water resources around the world, particularly over regions with insufficient in-situ measurements (e.g., Chen et al., 2019;Scanlon et al., 2018). Recently several studies (e.g., Hasan et al., 2018;Mehrnegar et al., 2020;Seyoum, 2018;Shamsudduha et al., 2017) have used GRACE observations to explore changes in TWS in the Nile River Basin (NRB).
Though, their conclusions are inconsistent because of different solutions for the GRACE observations used for investigative changes of TWS. Thus, it highlights the demand for evaluating the performance of each GRACE product in the NRB region.

Study area description
The Nile River Basin (NRB) is the world's longest river basin with a total flow length of over the year 2030, the population projected to rise to 700 million. According to (Abtew & Melesse, 2014), it raises concerns about sustainable and equitable management of the basin's resources.
Based on the hydro-climatic condition and recommendation of (Awange et al., 2014), we divided the NRB into four key areas of interest: (1) The Blue Nile Region (BNR) comprises a total area of 300,000km2. It mostly belongs to Ethiopian Highlands, including both the Upper and Lower Blue Nile region with Binda Rahad Region.
(2) Lake Victoria Region (LVR) has a total area of 258,000km2 representing Lake Kyoga, Albert, George, and Semilik basin. It is the headwaters of the White Nile.
(3) The Bahr-el-Ghazal Region (BER) that has a total area of 526,000km2. It represents the Congo-Nile River divide tributaries and the vast area of Sudanese plain with a low slope, which is the leading western tributary of the Nile. (4) Main Nile Region (MNR) has a total area of 258,000km2. Therefore, all NRB sub-basins (regions) exhibit more considerable spatial coverage that enables us to detect changes in TWS GRACE Satellite data. In this paper, the term region used to represents sub-basins of NRB interchangeably.

GRACE data
The GRACE data accessed from four official processing institutions. These are NASA Goddard Space Flight Center (GSFC), Center for Space Research at Texas State University (CSR), Geoforschungs Zentrum Potsdam (GFZ), and Jet Propulsion Laboratory (JPL). The processing institutions used different parameters and strategies, such as different background force models, data pre-processing approaches, local parameters, arc length, etc.
The TWS RL06 data of three SH solutions (CSR-SH, GFZ-SH, and JPL-SH) acquired from https://grace.jpl.nasa.Gov and http://geoid.colorado.edu/grace/. Numerous filtering methods used to remove north-south striping errors and correlated errors (Duan et al.,2009;Chambers & Bonin, 2012), most of which are empirical. In this study, we used the Level-3 GRACE data filtered by a de-striping technique to remove the correlated error (Swenson & Wahr, 2006). The SNR is already significantly improved (Göttl et al., 2018) in RL06 and posterior correction models (Fagiolini et al., 2015). The mass concentration (mascon) product is a recent advanced form of gravity field basis function. It is a much more rigorous approach comparing with the SH approach due to easy use for both geodesists and non-geodesists to implement the geophysical constraints without any postprocessing. In this study, we used three SH (CSR-SH, GFZ-SH, and JPL-SH) with their ensemble mean (Ens-SH) and three mascons (CSR-M, JPL-M, and GSFC-M) solutions (see Table 1) in a total of 7 solutions.

Ensemble modeling
The ensemble model could reduce the noise seen in the scattering of individual solutions.
According to (Sakumura et al., 2014), a weighted mean of the ensemble components will result in an over-all solution with a much higher level of reliability than a stand-alone model. Each GRACE solution is assumed to represent the truth plus a certain amount of error representing a relatively uniform scatter about this mean value. In this study, the ensemble model produced by composing the weighted mean of three SH GRACE solutions.

The amplitude of the GRACE TWS
The amplitude of the GRACE TWS signal was affected by different processing techniques. The smoothing filters and "de-striping" were applied to SH solutions by GRACE processing centers to reduce the noise level, which could reduce and smooth the amplitude of the water mass variations (Wouters et al., 2014). They also noted that the processing treatments and the truncation of the spherical harmonic degree suffer from leakage error. The leakage error impacts water mass that tends to be spatially spread and can "leak" toward nearby water mass. The global grid of timeinvariant scaling factors computed by (Landerer & Swenson, 2012) using the CLM4 global hydrological model could overcome such leakage to some extent. Details of these data (three SH and mascon solutions) used in this study are given in Table 1.

Rainfall data
Rainfall in the Nile River basin has been recorded since the beginning of the 20th century using the manual as well as automatic recording rain gauges. In this study, we used the rain gauge data collected from various sources. It includes respective countries' hydro-meteorological departments in the basin as well as from researchers in the region. Furthermore, to fill the missing gap and basin-wide representation, we had used Tropical Rainfall Measuring Mission (TRMM) 3B43 precipitation V7 data (1998 to 2016). TRMM satellite precipitation data is available and downloaded from the NASA website (https://trmm.gsfc.nasa.gov/). We used the TRMM-3B43 ( (see Table 1) Level 3 satellite-gauge (SG) combination data set to estimate the precipitation over NRB. The weighted areal mean technique is applied to determine the rainfall magnitude at the subbasin (regions) and basin level.

Groundwater storage
Groundwater storage ( ∆ ) changes are estimated in regional-scale (Richey et al., 2015;Ramillien et al., 2014) by computing total terrestrial water storage anomalies (∆ ) from one or more GRACE products and deducting changes in other terrestrial stores simulated by LSMs. In the NRB, it comprises changes in soil moisture (∆ ) and surface water (∆ ) storage change (Eq. (2)).
(2)), we used simulated soil moisture to represent ∆ and surface runoff, as a proxy for ∆ as presented in  after (Thomas et al.,2017); Mishra et al., 2016) that applied to derive from LSMs of NASA's Global Land Data Assimilation System (GLDAS). GLDAS is an uncoupled land surface modeling system that includes multiple global LSMs (Rodell et al., 2004), driven by the observed surface atmospheric fields (Sheffield et al., 2006). We apply monthly ∆ and surface runoff data at a spatial resolution of 1⁰

Satellite Altimetry of Major Lakes
The Water level estimation by satellite altimetry developed and optimized for open water bodies (Moore & Williams, 2014;Uebbing et al., 2015). In this study, NRB major lakes satellite altimetry-derived surface water observations used to validate GRACE products TWSA. The In this study, we used satellite altimetry-derived surface water observations to evaluate the performance of GRACE TWS ( Table 2). Details of these data are given in Table 1.

Soil moisture changes
Data from three different sources (GLDAS, WGHM, and ERA-Interim) merged using the Triple Collocation Analysis (TCA) (Gruber et al., 2017) after (Stoffelen, 1998) for reliable estimates of soil moisture changes over the NRB (see Figure 2). TCA offers an alternative method for estimating random error variances (Gruber et al., 2017) in the absence of ground reference data.
TCA is applied here to merge soil moisture outputs; with being the actual soil moisture variation, 1 , 2 , and 3 represents three soil moisture anomalies related to St with 1 , 2 and 3 being the coefficients that correspond to the errors of 1 , 2 , and 3 , respectively. The objective is to estimate error variances associated with 1 , 2 , and 3 to be used in the weighting process. On the one hand, TCA solves this by considering the errors of the products to be independent of each other. On the other hand, it arbitrarily assumes any of the products as a reference ( see Yilmaz et al., 2012) for more details regarding TCA' implementation.
Supposing the errors of the products are independent of each other and from the truth, and assuming a mutual linear relationship between the actual soil moisture and these estimates, the error variances ( 2 ) are calculated based on Eqs. 6-8: Selecting any of the products as the reference will not no impact the merged time series (Gruber et al., 2017). The error variances calculated from Eqs. (6) -(9) used in Eqs. (9) -(11) to estimate weights of each merged product where σ 1 2 , σ 2 2 , and σ 3 2 are error variances of S 1 , S 2 , and S 3 , respectively, with the corresponding weights of w1, w2, and w3.

Time series decomposition and trend analysis
The GRACE TWSA monthly time series (TS) data needs seasonal trend decomposition to remove seasonality from the original TWSA series (Andrew et al.,2017a) after filling missing data. The missing data can be filled by a simple temporal interpolation method (Andrew et al., 2017b) or by averaging the values of the two or more months either side. The harmonic analysis is a commonly implemented approach (Scanlon et al., 2016) and the seasonal trend decomposition using Loess (STL) exhibited similar decomposed results with harmonic analysis (Scanlon et al., 2018). In this study, the STLplus approach of LOESS (LOcal regrESSion) smoothing that is modified by (Ryan Hafen, 2016), which was initially introduced by (Cleveland et al., 1990), is used due to adaptability and robustness to decompose time series (see Eq. (13)).
Where it consists of inner and outer loops with a sequence of smoothing operator LOESS and generates three components from a time series; the original signal ( ) represented as the sum of a long-term part ( ), a seasonal cycle ( ), and the remaining sub-seasonal residuals ( ). The long-term component ( ) further divided into long-term linear trends and long-term nonlinear (inter-annual variability). The residuals reflect both sub-seasonal signal and noise (Humphrey et al., 2016); these high-frequency residuals anticipated being a combination of both the noise and a real signal that represents sub-seasonal water storage variability that is present in the GRACE data. Details of the STL decomposition approach presented in (Lu et al., 2003).
Thus, the TWS trends refer to the linear trends estimated from the long-term deseasonalized) after STL analysis. Moreover, the modified Mann Kendall Trend Test (M-K test) is used to identify patterns from the deseasonalized TWSA series. Yue et al. (2002) recommend the modified M-K test for a nonparametric and vigorous test for determining a trend in a time series. Detailed procedures presented on the respective reference given in appendix 2.

Uncertainty Estimates of GRACE solutions TWS
Uncertainties of all solutions estimated in this study are the Root mean square (RMS) of residuals of TWSA. The residuals are reminders after removing long-term trends, interannual, and seasonal signals ( i.e., ≥ 13 months) per the recommendation of (Scanlon et al., 2016). It may overestimate actual uncertainties since the residuals may contain sub-seasonal signals in accumulation to noise.
Leakage errors calculated using synthetic simulations as described in  and the shape of the river basin imitates to the placement of the mascons, associated with the spatial distribution of mass at sub-mascon spatial scales according to geophysical models, which were considered by leakage errors. These errors from the differences between original TWSA from the CLM4 LSM and the filtered TWSA applied with the scaling factors (Landerer & Swenson, 2012).

Results and Discussion
This section presents the data analysis results and provides a brief discussion about TWSA time series, trend (GWS, SWS, SMS), and uncertainty over NRB and sub-basin. In sequential order time series over NRB, GRACE TWSA trends and uncertainties over NRB discussed below.
Finally, the results of the spatiotemporal analysis of trends and uncertainty in the sub-catchments of NRB presented. better than other GRACE products (see details in Table 2), and thus it shows consistent with our study.

3.1.Time series over NRB
The overall seasonal cycle is captured very accurately and similarly for all GRACE TWSA except variations between the available solutions that strongly impact the hydrological information within the data. Sakumura et al. (2014) noted the solutions all varied from the mean 1-1.5cm level regardless of basin size over 156 global river basins. Thus, these differences among the solutions arise signal or extent of the regional signal and support the hypothesis that inconsistence among the solutions is due to random error rather than signal differences.

GRACE TWSA trends and uncertainties over Nile River Basin
GRACE TWSA time series data decomposed by using the STLplus techniques in terms of the long-term trend and seasonal TWSA signals ( Figure A1), and the residual presented in (Figure A1 (a-d)). The seasonal signals removed by STLplus deseasonalization performed in (Appendix 1 Figure A1d shows the capacity of STLplus to remove the seasonal variations from the original GRACE TWSA time series. According to the derived signals, the seasonal pattern of TWSA is similar to that of precipitation in NRB. It is lower in winter (dry season) and higher in summer (primary rainy season). Therefore, summer precipitation dominates the annual Rainfall in NRB (Onyutha & Willems, 2015), and abundant rainfall brings an increase in TWSA and water reserves in summer. The

Groundwater result
The result of TCA shows the average estimated soil moisture variation from GLDAS, WGHM, and ERA-Interim over the Nile regions Figure 7. Moreover, using the soil moisture and TWS changes, groundwater changes calculated based on (Eq. 1 & 2) are also presented in Figure 6 -7. Over different regions on NRB groundwater changes exhibit short-term and long-term trends.
Predominantly, over the Lake Victoria and Blue Nile Region, the negative trends can be seen between (2003-2006 and 2007-2009), followed by notable increases, probably due to the similar reasons clarified earlier studies (Awange et al., 2008;Seyoum, 2018). Pre-eminently, a negative groundwater trend is prevailing over the Main Nile River Region during the entire study period regardless of precipitation trends. After 2008 the rate of this decline is found to be larger, which shows substantial groundwater depletion over the region. The impact of the 2007 ENSO and reducing water controls in the Lake Victoria Region after 2006 likely caused groundwater increase over the Blue Nile Basin and Lake Victoria Region.
Though, this effect found to be tumbledown by 2008 follow-on in negative trends (with a higher rate for the Main Nile River Region) between 2008 and 2012. The groundwater to rise in both Lake Victoria and Blue Nile Region due to the increasing amount of precipitation after 2012.
Similarly, the high rainfall has the same impact on soil moisture variation between 2012 and 2019.
The spatial pattern maps generated (Figure 7) for three water storage compartments such as surface across MNR. This increase is due to the presence of significant surface Lakes such as Lake Tana, Lake Victoria, Lake Albert, and Lake Koga.
Similarly, a more significant share of LVR and BNR experienced increasing trend due to the presence of major natural lakes and reservoirs in the equatorial Lake Region and Lake Tana while in MNR region the decreasing trend (0.1 to 0.6cm/year) of SMS, SWS, and GWS due to intensive irrigation.

Main Nile Region (MNR)
The main Nile River (MNR) region is located in the northern part of NRB. It covers significant water resources such as Lake Nasser, Aswan high dam, and Nile delta in the Western Plateau within the Nubian Aquifer. The GRACE's TWS (Figure 8 & 9d) indicates a declining trend of stored water, while discrepancy among the GRACE solutions in terms of magnitude and pattern of the trend is very high. In the MNR region, the water lost consistently ( Figure A3d ) that may be due to two factors (Sultan et al., 2013). The first factor is that most of the water extracted from the Nubian Aquifer used for agricultural purposes such as the East Uweinat project that has seen heavy utilization of groundwater. In the East Uweinat project, the lands reclaimed amounted 1200 ha(1992), 4200 ha (2003), 75,000 ha by 2022, all of which will irrigate using groundwater as depicted by (Sultan et al., 2013). The second factor is the fact that the Uweinat-Aswan uplift prevents recharge of groundwater flowing from the South to the North.

Lake Victoria Region (LVR)
Lake Victoria Region (LVR) located in the southern part of NRB, and it is a source of the White Nile comprising Lake Victoria, Lake Albert, and Lake Tanganyika. In LVR (see on Figure A3b onward, it increases with a rate of 0.14cm/yr, and the water storage gain is 9.2cm in equivalent water height. The rising and declining trend show reasonable agreement with (Swenson & Wahr, 2009;Awange et al., 2008). The GRACE TWSA decrease between 2003-2006 is highly associated with the expansion of the Owen Falls dam/ Nalubale dam (Swenson & Wahr, 2009). On the contrary (Becker et al., 2010) reported an increase in precipitation from the end of 2005 due to the El Niño-Southern Oscillation (ENSO) effect explained as anthropogenic and climate variability contributions. The individual water balance terms retrieved from satellite altimetry by (Vanderkelen et al., 2018) indicate the increase of water level in the LVR from 2011 onwards.

Bahr-El-Ghazal Region (BER)
The Bahr El Ghazal Region (BER) represents rivers that originate from the Congo-Nile River divide in the Eastern part of NRB. Nearly all (i.e., about 96% ) the BER runoff and precipitation is evaporated or leaked into swamps (Conway & Hulme, 1993).  (Figure 9c), the overall trends statistically significant (p-value < 0.05) and are slightly rising (+0.24 to +0.13 cm/yr) except for the decline of JPL-SH. These discrepancies among the solutions exhibited higher over MNR (Figure 9d), which shows a slightly increasing for JPL-SH and the CSR-SH. On the contrary, for the rest of SH and mascon solutions, the decreasing trend noted. Similar trends are reported by (Shamsudduha et al., 2017)  The deviation from the ensemble model and variation of the trend among solutions are unsatisfactory to describe the relative quality of TWSA solutions of each processing center.
Henceforward, the signal-to-noise ratio (SNR) for each GRACE solution estimated as the ratio between the uncertainties and trend of the time series of TWSA.
Over the NRB, CSR-M performed best with the highest SNR and lowest uncertainty (see Table 2) that is followed by GSFC-M, while GFZ-SH showed the most inferior performance. Based on correlation and their Root mean square error (RMSE) with satellite Altimetry also CSR-M show the highest correlation (R2=0.87) and lowest (RMSE=0.52) followed by GSFC-M and JPL-M, respectively. In the NRB sub-basins (MNR, LVR, and BER) also mascon solutions performed very well while all SH performs the poorest.

Conclusion and Recommendation
This study compared and discussed different GRACE TWS solutions and their effects on the spatial-temporal variations of GWS, SMS, and SWS estimated in the NRB, which implicit that one should cautiously select suitable TWS solutions for future studies. The long-term TWS compartments (groundwater, soil moisture, and surface water storage) simulated using land surface models as well as multi-mission satellite products. Complimentary consistency between the CSR-M TWS and the satellite altimetry over the NRB that fully demonstrated the potential of GRACE data in monitoring changes in water (GWS, SWS, and SMS) dynamics over the NRB. Nevertheless, the extended groundwater level monitoring networks yet anticipated to convey more reliable groundwater information for validating the GRACE-based estimates furthermore. The key outcomes of the study summarized below.
1. GRACE derived TWS trends over the NRB and sub-basins evidenced the TWS time series deviations, uncertainties, trends discrepancies, and SNR variations among GRACE solutions of both mascons and SH. Consequently, CSR-M noted the first best performance; GSFC-M is the second-best, followed by JPL-M over NRB and at sub-basins levels. Thus, we recommend CSR-M solutions for detailed water resource studies to minimize uncertainties, especially on hydro-climatological extremes (droughts and floods).
2. Over the NRB, a significant long-term (2003-2017) negative groundwater and soil moisture trend found in the MNR (Sudan and Egypt), some part in LVR sub-basin (Tanzania & Kenya) demonstrating a potential depletion. Despite an increase in precipitation and TWS time series, the rate of decline seen to increase rapidly from 2008, thus indicating the possibility of human-induced decline ( e.g., for irrigation purposes).
Over the BNR, BER, and LVR sub-basins, negative trends found between January 2003 and September 2006.
3. More considerable spatiotemporal variations of TWS and precipitation recorded over LVR and a lesser degree over the BNR, which can explain more substantial water storage in these regions compared to the MNR (Egypt and Sudan). Contrary to the LVR and BNR, negative trends found for TWS variations over the MNR. In general, negative patterns are observed for water level variation in the MNR and BER regions with the highest (for the period 2003-2006) in contrast to the LVR and BNR.
Availability of data and materials: Data and materials will be made available on request.

Competing interests:
The authors declare no conflict of interest.

Funding:
We thank the National Natural Science Foundation of China (grant number 41974013 and 41574018) for financial support.
Author Contributions: Zemede M. Nigatu conducted the data analysis and wrote the paper; Dongming Fan and Wei You supervised, verified the methods and results, and revised the manuscript.