A Multi-Sourced Data Retrodiction of Remotely Sensed Terrestrial Water Storage Changes for West Africa

Abstract

Remotely sensed terrestrial water storage changes (TWSC) from the past Gravity Recovery and Climate Experiment (GRACE) mission cover a relatively short period (≈15 years). This short span presents challenges for long-term studies (e.g., drought assessment) in data-poor regions like West Africa (WA). Thus, we developed a Nonlinear Autoregressive model with eXogenous input (NARX) neural network to backcast GRACE-derived TWSC series to 1979 over WA. We trained the network to simulate TWSC based on its relationship with rainfall, evaporation, surface temperature, net-precipitation, soil moisture, and climate indices. The reconstructed TWSC series, upon validation, indicate high skill performance with a root-mean-square error (RMSE) of 11.83 mm/month and coefficient correlation of 0.89. The validation was performed considering only 15% of the available TWSC data not used to train the network. More so, we used the total water content changes (TWCC) synthesized from Noah driven global land data assimilation system in a simulation under the same condition as the GRACE data. The results based on this simulation show the feasibility of the NARX networks in hindcasting TWCC with RMSE of 8.06 mm/month and correlation coefficient of 0.88. The NARX network proved robust to adequately reconstruct GRACE-derived TWSC estimates back to 1979.

1. Introduction

A good understanding of regional hydrological conditions is vital to scientists, water resource managers, policy makers, among others, in helping to understand the state and flux of available freshwater. Such conditions are mostly evaluated with hydrological models, which enable the estimation and prediction of discharge and the various components of the terrestrial water storage. In West Africa (WA), the use of such models is nonetheless challenged by data scarcity [1]. Alternatively, global land surface models, such as, the Global Land Data Assimilation System (GLDAS) could be used to characterize hydrological conditions in the region [2]. However, these models perform poorly over the region due to the persistent issue of data deficiency [3]. The terrestrial water storage (TWS) inverted from the observations of the Gravity Recovery and Climate Experiment (GRACE) mission, therefore, is a very useful dataset which can be used to access changes in water availability to ensure a sustainable use [4]. However, the previous GRACE mission covered a relatively short period (April 2002 to July 2017), with 2 to 6 months latency, which makes it difficult to be used for long-term studies. Additionally, a data gap is envisaged between the past and the current (follow-on) GRACE missions. This, therefore, necessitates the need for an approach to forecast/backcast GRACE-derived TWS over the region.

As the variability of the water resources in WA is highly coupled to land-atmosphere-ocean interactions, both at the regional and global scales [5], it is reasonable to use hydro-climatic variables, such as rainfall, temperature, evaporation, and El Niño Southern Oscillation (ENSO) index, as retrieved from space-borne sensors, reanalysis datasets etc. to forecast/backcast TWS in the region. Using GRACE-derived TWS and in-situ river level records, Becker et al. [6] employed principal component analysis, which makes a linear and stationary assumption of time series to reconstruct GRACE-TWS from 1980 to 2008 over the Amazon basin. Similarly, de Linage et al. [7], developed a simple linear model between TWS over the Amazon and sea surface temperatures (SSTs) indices derived from Niño 4 and the Tropical North Atlantic Index (TNAI). Meanwhile, Forootan et al. [5], developed an autoregressive (ARX) model, which utilized independent component analysis (ICA) for dimension reduction to predict TWS over WA using rainfall data from the Tropical Rainfall Measuring Mission (TRMM), and SSTs of the three major oceanic basins (Pacific, Indian and Atlantic). However, this ICA/ARX method assumes a stationary state of TWS over the region and the prediction accuracy reduces after a forecast period of two years. Furthermore, Yin et al. [8] used the water balance method to extend the TWS series back to 1980 using multi-source datasets as inputs in the terrestrial water balance equation.

All the examples herein given require the use of mathematical models. It must be noted that, since the formulation of these models are based on empirical data, it essentially makes this a machine learning/pattern recognition problem [9]. Mukhopadhyay [10] indicated that no mathematical model could effectively characterize hydrological phenomena. Thus, algorithms which make no prior assumption(s) of time series with adaptive capabilities offer a good alternative for predicting TWS over a region. Artificial Neural Networks (ANNs), which are effective machine learning tools, offer such algorithms that are, self-learning, self-adapting, and self-organizing and capable of predicting hydrological variables with high efficiency [11]. Long et al. [12] demonstrated this by training an ANN, which was used to extend GRACE-derived TWS series back to 1982 over the Yun–Gui Plateau and its sub-regions, China.

The ANN machine learning algorithm is exceptionally well suited for modeling input–output relationships, especially, in the absence of optimally calibrated physically-based models [12]. ANN has been extensively used in forecasting hydro-climatic parameters such as, stream flows, groundwater, rainfall, and droughts, with reasonably great performances (see References [13,14,15]). The major strength of ANN lies in the ability to efficiently learn the causal relationships within a nonlinear dynamic system, without a priori assumption(s) of the physical processes in the system [16]. In the context of GRACE-derived TWS, there is a growing interest in the use of ANNs, although the applications are relatively few. Yirdaw-Zeleke [17] applied a neural network to downscale GRACE-derived TWS into local groundwater storage over the Assiniboine Delta Aquifer, Canada. Likewise, Miro and Famiglietti [18] downscaled the GRACE-TWS to high-resolution groundwater storage in California’s Central Valley. Sun [19] employed it to predict groundwater level using GRACE data over the mid-west regions of the USA. Moreover, ANNs have been used to extend GRACE data to the early 1980s over mainland Australia [20] and karst plateau regions in southwest China [12]. Notable are also the contributions of Zhang et al. [21], who extended TWS series over Yangtze basin, and Mukherjee and Ramachandran [22], who used ANN (among other algorithms) to predict groundwater variations in India. Chen et al. [23], monitored the hydrological extremes (droughts and floods) Liao River Basin, Northeast China, using extended GRACE data based on ANN.

To the best of our knowledge, this current study represents the first attempt to apply ANNs in the context of TWS changes (TWSC), derived from GRACE over any region in WA. Additionally, we make use of a broad data spectrum, which replicates dynamics in the water and energy cycles that impact TWSC to train the network. In this study, we present a nonlinear Autoregressive Neural Network with eXogenous Inputs (NARX) which was used to backcast GRACE-derived TWS over WA. The exogenous data inputs employed included: rainfall, evapotranspiration, land surface air temperature, precipitation minus evapotranspiration (from both the atmospheric and land perspective), soil moisture as well as global and regional circulation indices. The NARX was used to learn the physical relationships between GRACE estimates and the afore-mentioned fields over the period 2003 to 2013. The trained network was then used to retrodict TWSC estimates over WA from 2013 to 1979 (34-years period). To this end, for the first time, we consider the accuracy of backcasted TWSC from synthetic GRACE-TWSC to address the question of whether extended TWSC series can be used to address the past seasonal variations in TWSC over WA.

The rest of the study is outlined as follows: the datasets and the method used are presented in Section 2; while the Section 3 describes the results followed by their discussions in Section 4. Summary of the study is presented in Section 5.

2. Materials and Methods

This section therefore briefly presents the specific data employed to backcast the GRACE-derived TWS as well as the reasons behind their usage over the study area in West Africa (Figure 1). All datasets and methods used in this work are described in detail in Section 2.2 and Section 2.3. A brief description of the ANN deployed for the TWS forecasting is also given here.

2.1. Study Area

West Africa covers an area of approximately of equivalent to 20% of the total area of Africa lying between longitudes 18° W and 16° E and latitudes 3° N and 28° N (Figure 1). The region is bounded in the west and south by the Atlantic Ocean and in the north by the Sahara Desert. The region’s most important highlands include the Guinea Highlands (including the Fouta Djallon, Guinea), which is the source of many rivers (e.g., Niger, Senegal, Gambia); encompassing parts of Guinea, Sierra Leone, northern Liberia, as well as northwestern Côte d’Ivoire; the Jos Plateau in Nigeria; Mandara Mountains, the source of the Benue river in southeastern Nigeria. The most significant river is the Niger, the fifth largest river in the world and second in Africa, which drains an area of 2 million km². Other important rivers in the region are the Gambia, Senegal, Comoe, and Volta.

The region is partitioned into sub-climatic zones using the K-means clustering algorithm. As suggested by Andam-Akorful et al. [24], this reclassification of the climatic zones considering their annual rainfall amounts to maintain consistency with the general African sub-climatic classification suggested in by [25]. Thus, the climatic zones are classified as hyper-arid (HA) in the north, arid, semi-arid (SA), dry sub-humid (DSH) and moist sub-humid (MSH) in the south (see Figure 1) the hyper-arid sub-region coincides with the desert area; arid, Sahel; semi-arid, Sudanian; Guinean, the moist and dry sub-humid areas [24]. West African rainfall is highly variable leading to the region experiencing recurrent droughts and floods.

2.2. Datasets

2.2.1. ARCv2 Rainfall

ARCv2 is a revised version of Climate Prediction Center’s (CPC) African Rainfall Climatology (ARC), which provides Africa specific precipitation from 1983 onwards. This product generates fields based on a subset of source data used in the Rainfall Estimate version 2 (RFE2) algorithm at a spatial resolution of 0.1° [26]. ARCv2 input dataset consist of input data from a 3-hourly geostationary infrared (IR) sensor centered over Africa by the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) and the Global Telecommunication System (GTS) gauge observations. Daily data retrieved from [27] were aggregated to monthly fields at a spatial resolution of 0.25°.

The product was chosen due to its reasonable performance over WA compared to other contemporaneous datasets, such as TRMM, PERSIANN, TAMSAT, and CMORPH [28]. Additionally, it is available from 1983, providing enough data for long-term studies. Note that, the ARCv2 rainfall fields are, however, not available from January 1979 to December 1982. Thus, in Appendix A we used a dedicated network to extend its records to 1979 using rainfall estimates for the Noah driven GLDAS version 2 (V2) land surface model [2], as well as, ENSO and Atlantic Niño indices. The predicted rainfall series showed a strong agreement with the original series over their common overlapping periods (see

Table A1. Area-averaged statistics for backcasted ARCv2 rainfall estimates over river basins and sub-climatic zones (see Figure 1 for the regions).

Region		ME (mm/month)	RMSE (mm/month)	NSE	R2
Basin	Niger	2.76	13.10	0.90	0.92
	Senegal	1.56	16.50	0.82	0.83
	Volta	5.36	23.01	0.83	0.84
Climate zones	Humid	−2.64	30.90	0.82	0.88
	Dry sub-humid	4.65	20.66	0.90	0.92
	Sudanian	4.81	18.47	0.89	0.83
	Sahelian	3.76	13.90	0.76	0.84
	Sahara	1.13	4.92	0.26	0.92
WA		2.21	11.00	0.89	0.92

and Figure A1).

2.2.2. ERA-Interim Data System

In order to model response of TWS to evapotranspiration and temperature changes, instantaneous moisture flux (equivalent to evapotranspiration [29]) and land surface air temperature (LSAT) fields were retrieved from ERA-Interim data system. Additionally, precipitation ($P$) minus evaporation ($E$), also known as net-precipitation, derived from the atmospheric water budget [24] was used to further characterize TWS response to atmospheric moisture conditions over the region, while soil moisture, also obtained from ERA-Interim was used to model land surface moisture conditions. The datasets were all obtained at a gridded spatial resolution of 0.25° and a monthly temporal interval from [30]. The motivation to predominantly select ERA-Interim reanalysis data for this study is due to its long-term consistency [31].

2.2.3. GLDAS-Noah Version 2

Additionally, total soil moisture (the sum of all layers), canopy and snow storages from GLDAS-Noah [2] in its Version 2 (V2) was also used to validate the result of the ANN predicted TWS, we named it TWC (total water content). GLDAS-2 datasets consist of 0.25°-by-0.25° gridded data, covering the period from 1948 to 2010. The temporal resolution for the GLDAS products is 3-hourly, where monthly products are generated through temporal averaging of the 3-hourly products. The data was retrieved from [32].

2.2.4. Climate Indices

The Bivariate ENSO Time series (BEST, [33]), which is derived from Niño 3.4 and Southern Oscillation Index (SOI) was used to model the response of WA to the ENSO phenomenon. The BEST index combines the atmospheric (SOI) and oceanic (Niño 3.4) components of the ENSO processes into a single field, and thereby, provides a more realistic characterization of the phenomenon. Furthermore, the Atlantic Equatorial Mode, also known as the Atlantic Niño [34], which modulates inter-annual rainfall over the region through the West African Monsoon (WAM) winds [35] was also used to simulate the regional oceanic-atmospheric circulation processes that drive precipitation patterns from the Gulf of Guinea. Finally, the time series of global temperature anomalies were also used to simulate the effect of global warming on TWS variability over the region as the droughts in the late 20th century have been linked to rising temperatures globally.

2.2.5. GRACE Fields

GRACE-derived TWS were computed at a spatial resolution of 1°-by-1° using level 3 products [36], from the three major processing centers—the Center for Space Research (CSR), Jet Propulsion Laboratory (JPL), and the GeoForschungsZentrum (GFZ). The monthly fields were scaled [36], and missing months in the GRACE series were obtained through a bilinear interpolation algorithm. GRACE solutions from the different centers (CSR, GFZ, JPL) present differences due to the unique processing methodologies of respective agencies [37]. In order to obtain of optimal estimates, Sakumura et al. [38] suggested the use of the simple average of the solutions from the different centers. However, since the respective products perform differently in different parts of the world, an uncertainty-based weighting of the ensemble product will offer a more equitable solution [39]. In this study, a different approach which utilizes uncertainty-based weighting is proposed to build an ensemble model. To this end, the Three-Cornered-Hat (TCH) was employed to evaluate the uncertainties of estimates from the different centers in the absence of any reference data. Details are provided in [39] and MATLAB^® (The MathWorks, Inc., Natick, MA, USA) scripts are available under request to the authors. The uncertainties were used as weights in the realization of the ensemble model using the expression:

$${\mathrm{TWS}}_{\mathrm{ENM}}\left(t\right)={{\displaystyle \sum }}_1^3{w}_i\left(t\right){\mathrm{TWS}}_i\left(t\right),$$

(1)

where $w_i\left(t\right)$ is the time-dependent normalized weight based on the uncertainties $\delta$, as obtained from the TCH, expressed as:

$$w_i=\frac{\frac{1}{{\delta }_{{\mathrm{TWS}}_i}^2}}{{{\displaystyle \sum }}_1^3\frac{1}{{\delta }_{{\mathrm{TWS}}_i}^2}}.$$

(2)

Additionally, with the knowledge in the uncertainties within the respective series, the signal-to-noise-ratio (SNR) of each product is computed as a ratio of the RMS to the noise estimate. The purpose of this metric is to assess the amplitude of the TWS field concerning the uncertainties in order to gain a better understanding of the amount of “contamination” within GRACE-derived TWS.

We furthered expressed TWS fields in terms of TWS changes (TWSC, i.e., fluxes) using numeric differentiation using the central derivative operator as:

$$\mathrm{TWSC}\left(t\right)=\frac{\mathrm{TWS}\left(t+1\right)-\mathrm{TWS}\left(t-1\right)}{2},$$

(3)

where $t$ is the month. This filter is important in order to compare the TWS with other fluxes, or, in our case, use other fluxes (e.g., rainfall) to backcast TWSC. We also will backcast TWS for comparison with previous works [21].

2.3. Methods

The methodology presented in the next subsections are summarized in Figure 2. First, we constructed a robust time-series for TWS (TWS_ENM), considering the products computed from the solutions provided by CSR, GFZ, and JPL (see Section 2.2.5). Second, the artificial neural network was considered in a two-fold experiment. Simulation of TWS (TWC) from our preferred hydrological model is used in a “closed-loop” simulation to assess the quality of the extended TWS series from 2003 to 1979 (Section 2.3.1). Third, the outputs from the artificial neural network are validated using the metrics presented in Section 2.3.2.

Due to the lack of long records for the purposes of evaluation, the backasting approach was further validated by designing a similar system to reconstruct TWC from the Noah V2 LSM, which is endowed with long-term data series. This network was trained to predict TWC from 2001 to 2010, after which it was used to backcast TWS from 2010 to 1979. Records from 1979 to 2000 from the reconstructed and the original datasets was then used to validate results from the ANN (Figure 2).

2.3.1. Reconstruction of GRACE-Derived TWS/TWSC

For this study, the Neural Network with eXogenous Inputs (NARX), which is suitable for modeling the nonlinear system was used in the design of the network. It is a dynamic neural network and contains recurrent feedback from several layers of the network to the input layer [40]. A NARX network predicts a signal by regressing the previous values of the output signal and previous values of an independent (exogenous) input signal. Mathematically, a NARX network is expressed as [40]:

$$y\left(n\right)=f\left[y\left(n-1\right),y\left(n-2\right),\dots ,y\left(n-n_y\right),u\left(n-1\right),u\left(n-2\right),\dots ,u\left(n-n_u\right)\right],$$

(4)

where $y\left(n\right)$ is the signal to be predicted, $u\left(n\right)$ are the independent exogenous input signal at a time discreet time step $n$, $n_u$ and $n_y$ is the input and the output delay, and $f(·)$ is the transfer function (normally unknown and can be approximated). A NARX architecture is shown in Figure 3 with one time series as input and another as output. In Figure 3, $n$ is the discreet time step, $h$ is the number of hidden layers, $w_i$ is the weight of input layer, $w_0$ is the weight of the output layer, and $v$ is the parameter of the $i$-th layer.

Our backcasting model employed 16 hidden layers and calibrated with a Bayesian regularization back propagation learning algorithm. Furthermore, the sigmoid transfer function was used. The exogenous signals consisted of the several fields, including: rainfall, evapotranspiration, land surface temperature, soil moisture, net-precipitation (P − E) from the atmospheric and terrestrial water budget approaches, and climate indices (BEST and AEM). All the gridded datasets, apart from the GRACE were resampled to a spatial resolution of 0.25° from 1979 to 2013 (34-year period).

Upon readying all needed exogenous variables, 2003 to 2013 was selected as the training period. The trained network was then used to reconstruct TWS from 1979 to 2013. In the designing the network, 70% of the data was used for training, whereas 15% for validation and to stop training before overfitting; and 15% for testing (used as independent data). Additionally, GRACE-derived TWS solutions covering the year 2002 were used to validate the results obtained from the ANN. The same experiment was performed for TWSC, that is, using flux (TWSC) rather than storage (TWS). The resulting solutions from the ANN were smoothed with a discrete cosine transform based smoothing algorithm [41]. We used the Neural Network Toolbox [42] of MATLAB^® and the script is available under request to the authors.

To further validate the GRACE estimates in the face of limited data, a similar neural network was created for the GLDAS-Noah simulated total water content (TWC, a combination of all soil moisture layers and vegetation canopy water) data, which is endowed with a long-term series. This network was trained with similar exogenous variables as that of the GRACE-derived TWS (Figure 2) over a 10 year period, that is, 2000 to 2010. The network was then used to predict Noah-derived TWC estimates from 1979 to 2010, and the results compared to the original dataset over the period 1979 to 1999. Like for GRACE data (TWS and TWSC), we performed the same training however for TWCC (i.e., flux after applying the temporal derivative for TWC). This was necessary in order to validate the retrodicted results over the time period of 1979–2000 (for GLDAS-Noah).

2.3.2. Validation Metrics

All backcasted series (ARCv2, GRACE-derived TWS, and Noah) were all evaluated using four validation metrics, which included, the coefficient of determination ($R^2$); mean error (ME); root mean square error (RMSE); and the Nash-Sutcliffe efficiency coefficient (NSE).

The coefficient of determination allows the quantification of the common information content between observed and backcasted series and is given by:

$$R^2=1-\frac{{{\displaystyle \sum }}_1^n{\left(f_i-\overline{y}\right)}^2}{{{\displaystyle \sum }}_1^n{\left(y_i-\overline{y}\right)}^2},$$

(5)

where $f_i$ is the backcasted value at an epoch, $I$ and $y$ is the original value. The range of $R^2$ lies between 0 (no correlation) and 1 (perfect fit), which describes how much of the observed dispersion is explained by the backcasting. The drawback of $R^2$ when considered alone is that, a model with large positive/negative bias may still show a good correlation and thus, require the additional consideration of other metrics such as ME and RMSE.

Alternatively, the NSE is given as:

$$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_1^n{\left(y_i-f_i\right)}^2}{{{\displaystyle \sum }}_1^n{\left(y_i-\overline{y}\right)}^2}.$$

(6)

NSE ranges between $-\infty$ and 1, 1 signifying a perfect fit (the backcasted and original estimates agree in mean, amplitude, and phase), whereas an efficiency lower than zero indicates that the mean value of the observed time series would have been a better predictor than the model (the ANN in this case).

Like $R^2$, NSE is not very sensitive to model over- or underprediction, especially during dry epochs [43]. RMSE is obtained as:

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_1^n{\left(f_i-y_i\right)}^2}{n}}.$$

(7)

Finally, ME is given as

$$\mathrm{ME}=\frac{1}{n}{{\displaystyle \sum }}_1^n(f_i-y_i).$$

(8)

In Equation (8), if dividing the residuals by the observed field $y_i$ would provide the mean bias error (MBE).

3. Results

3.1. The Ensemble GRACE-TWS Fields

Figure 4 illustrates the spatial distribution of the uncertainties and SNRs based on TCH method. It is evident that CSR estimates present lower uncertainties relative to GFZ and JPL (compare Figure 4a–c, respectively).

The spatial distribution of the uncertainties shows high estimates along the coastlines, especially around North America and Russia, whereas low estimates are observed in arid regions, for example, the Sahara and Gobi deserts (Figure 4a–c).

On the other hand, high SNRs are observed in places with high amplitudes in TWS, such as, the Amazon basin (Figure 4d–f). Again, CSR (Figure 4d) shows the best performance as expressed by SNR values.

Table 1. Distribution of uncertainties and SNRs for CSR, GFZ, and JPL TWS. The percent shows the percentage of cells lying between (Lower Boundary—Upper Boundary) for each product.

Product	Uncertainties Intervals (mm)
	0–20	20–40	40–60	60–80	80–100	100–120
CSR	91.5%	7.4%	0.8%	0.2%	0.1%	0.0%
GFZ	63.6%	29.7%	4.9%	1.3%	0.4%	0.1%
JPL	77.9%	18.7%	2.6%	0.6%	0.1%	0.1%
Product	SNR Intervals (Unitless)
	0–20	20–40	40–60	60–80	80–100	100–120
CSR	95.5%	3.7%	0.5%	0.2%	0.1%	0.0%
GFZ	99.9%	0.1%	0.0%	0.0%	0.0%	0.0%
JPL	99.5%	0.5%	0.0%	0.0%	0.0%	0.0%

presents the distribution of the magnitudes of uncertainties and SNRs, expressed in percentages over landed areas globally with 91.5% of its estimates falling within 0–20 mm of the noise. Similarly, CSR presents slightly higher SNRs as opposed to the other two products.

The comparisons between ensemble models, constructed with simple averaging, that is Equation (1) without weights, and uncertainty-based weighting, that is Equation (1) considering the weights given by Equation (2), are given in

Table 2. Distribution of uncertainties and SNRs for the simple ensemble average (ENM_AVE) and Three-Cornered-Hat (TCH)-weighted ensembles (ENM_TCH). The percent shows the percentage of cells lying between (Lower Boundary–Upper Boundary) for each product.

Product	Uncertainties Intervals (mm)
	0–10	10–20	20–30	30–40	50–60	60–70
ENMAVE	87.1%	11.2%	1.2%	0.3%	0.1%	0.1%
ENMTCH	97.1%	2.6%	0.2%	0.1%	0.0%	0.0%
Product	SNR Intervals (Unitless)
	0–20	20–40	40–60	60–80	80–100	100–120
ENMAVE	85.7%	11.1%	1.6%	0.9%	0.3%	0.1%
ENMTCH	64.1%	27.4%	5.3%	1.8%	0.9%	0.5%

.

About 97% of the cells distribution, present uncertainties lower than 10 mm in the weighted model, ENM_TCH, as opposed to 87% in ENM_AVE (Table 2). Although both products have higher SNRs than the original series from the processing centers, the weighted product (ENM_TCH) indicate higher SNRs compared to the averaged one (ENM_AVE). Whereas approximately, 86% of the cells in ENM_AVE have SNRs lower than 20, only about 64% of those of ENM_TCH fall within the same range, while significant 27% of the grids have their SNRs between the 20–40 range, in contrast to the 11% of ENM_AVE.

The multiple comparison procedure (MCP) [44,45] was also applied to evaluate the statistical differences between the original products and the two ensemble models, as well as, a statically based ranking of the means of the different datasets. The results of the MCP confirm that, the mean uncertainties in the combined models are lower than the original series (Figure 5). Again, the TCH weighted model outperforms the averaged one both in terms of uncertainties and SNRs, as shown in Figure 5. (Note, low values of mean uncertainties better the model, conversely, higher values of SNRs imply better the model.)

The procedure ranks the original products in the following descending order: CSR, JPL, and GFZ for both, uncertainties and SNRs. It must be noted here that, this is not to be used as a basis of choosing one product over the other, as many studies, including our own, indicate that products from the different processing centers essentially yield similar results. However, the use of a combined product, especially weighted by uncertainties is strongly recommended as shown the current study, in agreement with Sakumura et al. [38].

3.2. Confirmatory Study Based on Noah-Derived TWCC Prediction

In Figure 6, the results of the confirmatory validation by comparing ANN predicted Noah TWC and the observed estimates are presented.

Generally, the more humid south (apart from the very humid southwest) exhibited negative biases (Figure 6a), while the mid-latitude areas presented moderate to high overestimations. The hyper-arid areas in the north presented low RMSEs (Figure 6b), mostly less than 5 mm. The RMSEs for most other areas generally ranged between 10 to 20 mm, while relatively large values were found over wetter regions in south WA, as well as parts of the eastern section.

Correlations between the two datasets (predicted versus the actual) were high in the relatively humid south ($R^2\ge 0.7$) in contrast to the low correlation coefficients in the dry north (Figure 6c). At the high latitudes, some areas even presented negative correlations, indicating a likely weak performance of the network in the desert areas. The weak performance in this area is however of negligible effect to analysis performed in this study as most of the TWS variance occur from the Sahelian area down to the coastal south. Similarly, NSE, as shown in Figure 6d indicate high performance in the south relative to the desert areas in the north.

Table 3. Area-averaged statistics for backcasted Noah-derived TWCC over river basins and sub-climatic zones (see Figure 1 for the regions).

Region		ME (mm/month)	RMSE (mm/month)	NSE	R2
Basin	Niger	0.15	6.21	0.90	0.92
	Senegal	−0.96	7.70	0.83	0.83
	Volta	−1.58	10.40	0.88	0.90
Climate zones	Humid	1.24	10.90	0.89	0.92
	Dry sub-humid	1.18	8.42	0.93	0.94
	Sudanian	−1.04	12.46	0.84	0.86
	Sahelian	−0.95	6.78	0.74	0.76
	Sahara	0.06	1.09	0.51	0.57
WA		0.13	5.01	0.91	0.93

indicates high skill scores for all spatial averages. Thus, it can be concluded that, the ANN estimates are of high quality and can be used to adequately characterize spatial and temporal patterns in freshwater reserves over WA. (This is presented in the next, Section 3.3.)

3.3. Evaluation of Backcasted GRACE-Derived TWSC

The performance skills of the GRACE-derived TWSC are illustrated in Figure 7.

The MEs (Figure 7a) range from −20 to 40 mm/month across WA, although values range from approximately −10 to 10 mm/month in most areas. The RMSEs (Figure 7b) are lowest over the Sahara (≤15 mm/month), while the areas south of the Sahara generally present values between 20 and 50 mm/month. The highest RMSEs were found around the Niger Delta, with values around 100 mm/month.

The coefficients of determination (Figure 7c) indicate low correlations in the Sahara (R² ≤ 0.5), whereas the values in the Sahel and below are greater than or equal to 0.6. The spatial patterns of NSEs are like those of the R², with worst efficiencies being observed in over the Sahara, and values greater than or equal to 0.5 occurring from the Sahelian areas to the more humid south (see

Table 4. Area-averaged statistics for backcasted GRACE-derived TWSC over river basins and sub-climatic zones (see Figure 1 for the regions).

Region		ME (mm/month)	RMSE (mm/month)	NSE	R2
Basin	Niger	0.08	8.26	0.91	0.91
	Senegal	−0.50	10.68	0.78	0.78
	Volta	−0.01	13.23	0.89	0.89
Climate zones	Humid	2.40	14.83	0.93	0.84
	Dry sub-humid	0.47	14.17	0.91	0.91
	Sudanian	−0.22	13.00	0.86	0.78
	Sahelian	−0.67	5.16	0.81	0.89
	Sahara	−0.66	1.93	0.26	0.91
WA		0.05	6.98	0.91	0.91

and Figure 1).

All major basins and sub-climatic areas as presented in Table 4 show high skill score in all four metrics, apart from the Sahara’s NSE score (0.26). The highest RMSE was observed in the humid area, whereas the lowest was expectedly recorded for the Sahara.

Figure 8a is a plot of the WA spatially-averaged estimates of TWSC_ANN, TWSC_obs, and TWCC_Noah.

The strong agreement between the three series is evident. It is worth mentioning here that, TWCC_Noah was not considered as an exogenous data input during the training of the ANN used to retrodict GRACE-TWSC. Its strong agreement with the reconstructed GRACE series therefore, indicates robustness of the predicted estimates over the region.

Similarly to TWSC, Figure 8b presents spatially-averaged time series of GRACE-derived TWS estimates from the ANN (TWS_ANN) as well as the observed (TWS_obs) and Noah (TWC_Noah) estimates. We performed this in order to compare with previous results [21,23]. The three series show a strong agreement over their respective overlapping periods. TWS_ANN presents a high correlation with TWS_obs with an R² of 0.90 and an RMSE of 12.92 mm. The TWS_ANN signal amplitudes between 2003 and 2007 are slightly lower than those of TWS_obs, however, the converse is observed from 2008 onwards. With respect to TWS_Noah, an R² of 0.89 was obtained, indicating a high correlation. The Noah estimates tend to present a wet bias throughout the period of data availability (1979 to 2010). This seems to be because model-based estimates of TWS are less variable than GRACE-derived TWSC primarily due to the absence of underground water storage as well as surface water reservoir storage.

In Figure 8b it seems that inter-annual variations exist for the stored water through the decades. Thus, in

Table 5. Linear trends estimated from retrodicted GRACE and Noah series covering different time spans.

Product	Span	Linear Trend (mm/year)	p-Value
TWSobs	2002–2013	3.64 ± 1.20	3.8 × 10−3
TWSANN	2002–2013	3.44 ± 1.06	1.4 × 10−3
TWSANN	1979–2013	1.04 ± 0.21	6.5 × 10−7
TWCNoah	1979–2010	0.76 ± 0.19	7.5 × 10−7
TWSANN	1979–2006	0.67 ± 0.28	1.8 × 10−2
TWSANN	2007–2013	4.36 ± 2.38	7.1 × 10−2

it is presented the long-term trends for the series considering different time periods.

Generally, the backcasted TWS over the region presented an upward linear trend of 1.04 ± 0.21 mm/year from 1979 to 2013 (34-year period), whereas TWC_Noah indicated a positive trend of 0.76 ± 0.19 mm/year (Table 5). The linear trend from 2007 to 2013 is especially strong with a value of 4.36 ± 2.40 mm/year, while 1979 to 2006 presented a value of 0.67 ± 0.28 mm/year. All statistical measures were computed at a 95% confidence level.

4. Discussion

One of the main goal of this experiment was to attempt to assess the quality of backcasted GRACE-derived TWS over WA. We considered a “closed-loop” simulation by using synthetic TWSC (i.e., TWCC) predicted by a land surface model, GLDAS-Noah V2, since it covers the entire period of study. A period equivalent to the GRACE observation were used to generate an artificial neural network-based learning machine (ANN NARX) and then retrodict the TWCC till 1979. The network was trained to simulate Noah TWCC signals from 2003 to 2010 over WA, based on their physical nonlinear relationships with seven hydro-climatic variables (rainfall; evapotranspiration; land surface air temperature; net-precipitation; soil moisture; ENSO index; and global temperature anomaly). This was necessary in order to validate the ANN NARX. Overall, the results of this network presented similar spatial accuracies as that of GRACE, whereas the spatially-averaged TWCC presented ME, RMSE, NSE and $R^2$ of 0.13 mm/month, 5.01 mm/month, 0.91, and 0.93, respectively (Table 3). Furthermore, considering the main river basins over WA and the climate zones (cf., Figure 1), the errors are within GRACE-derived TWS [38,39]. There is some overestimation over moist sub-humid areas (Figure 6a) although NSE values shows reasonable performance. Yet, the ANN method seems to be robust to adequately reconstruct GRACE-derived TWSC estimates over WA as shown by the simulation.

Following this, the second goal was to retrodict the actual GRACE-derived TWSC to 1979 using a similar network. Firstly, an ensemble model of GRACE-derived TWS was created using the data from the three primary processing centers (CSR, JPL, and GFZ). In this regard, we applied the TCH method to independently estimate uncertainties in GRACE-derived TWS. The overall uncertainties estimated for the individual processing centers show that CSR delivered the lowest values. This is somewhat in agreement with previous studies [38,39], which found CSR providing the lowest noises among three processing centers. The computed uncertainties were subsequently used as weights in creating an ensemble TWS model from the original series, as well as, another combined model obtained through simple averaging. The weighted ensemble product, when compared to the original and the averaged datasets presented lower noise estimates and higher SNRs (Table 1 and Table 2), indicating an improvement in the TWS solution. This improvement in the ensemble series using weights contrasts with the results provided in Ref. [38].

Following this, the improved product was used in the reconstruction of GRACE TWSC (TWS was converted to TWSC by means of Equation (3)) retrodicted from 2013 to 1979 (34-years period) using an ANN NARX. Figure 8a depicts the averaged series for WA, which shows an overall agreement with observed TWSC with ME, RMSE, NSE, and $R^2$ of 0.05 mm/month, 6.98 mm/month, 0.91, and 0.91, respectively. Despite the differences in the length of the comparisons, the results of a “closed-loop” simulation (Noah-TWCC) and those for actual data (GRACE-TWSC) show consistency for the entire region and its sub-domains (compare Table 3 and Table 5). Furthermore, the use of TWSC shows relatively low NSE values in comparison with previous studies that adopted TWS [12,21,23], known that there are differences in the datasets, time span, study region, etc.

To support the above discussion that TWSC is a better candidate to be predicted than TWS an ANN NARX was trained to predict GRACE-TWS with the same hydro-climatic variables. Similarly, the ANN NARX predicted GRACE-derived TWS series (TWS_ANN), showed strong agreement with the original GRACE data, presenting similar trends over their overlapping periods (Figure 8b). The results, upon validation showed high skill scores, mostly in areas south of the Sahara Desert, but low performance in the desert due to the very low amplitudes in signals. The coefficient of determination and NSE coefficient values obtained for areas south of the Sahara were mostly greater or equal to 0.6 and 0.5, respectively, while RMSEs ranging from 20 mm to 50 mm were predominant. Spatially-averaged series for major basins and sub-climatic zones, as well as, the whole of WA itself presented median RMSE, NSE and $R^2$ of 11.83 mm, 0.76 and 0.89, respectively. This is somewhat in agreement with the findings presented in [12], thought different scenarios. Interesting, Figure 8a shows a better agreement between GRACE-TWSC with Noah-TWCC in comparison with the results depicted in Figure 8b, which shows GRACE-TWS and Noah-TWC. That means TWSC over WA mainly reflects fluctuations in soil moisture while TWS still have some memory effect, mainly in the groundwater compartments.

As mentioned in the Section 1 (Introduction), so far no one appear to have assessed ANN considering a “closed-loop” simulation. Thus, the importance of results using such algorithm lies both to validate the algorithm in the aim that the errors lie within an acceptable range (e.g., GRACE errors) and to provide long-term series to enable long term studies (e.g., droughts). However, we have not considered and investigated the contribution of desiccation of Lake Chad [46] on the hindcasted TWSC. Furthermore, the regularization of Lake Volta due to the water impoundment must be considered. Sensitivity of the output in ANN to the set of the parameters has not been considered. We just used seven hydro-meteorological variables that could modulate TWSC over WA, however, it is recommended to select inputs one-by-one or several at a time to assess the prediction performance to a set of inputs. Nonetheless, our finds support the use of ANN to extend TWSC series dedicated to long-term studies over West Africa, a data-poor region.

5. Summary

Although GRACE currently offers the most viable option of obtaining reliable estimates of TWS at regional and global scales, its time series is relatively short. For example, it cannot be used to infer the long-term changes in water availability over West Africa since the 1970s. Thus, a data-driven machine learning approach, which involved the use of NARX network was adopted to retrodict the GRACE series to 1979. The network was trained to learn the complex nonlinear relationships between GRACE-derived TWSC and the following hydrological variables: rainfall; temperature; evapotranspiration; net-precipitation; soil moisture, ENSO and Atlantic Niño (Niña) indices; and global temperature anomalies over the period 2003 to 2013. The trained network was subsequently used to backcast TWSC estimates from 2013 to 1979 covering a period of 34 years. Due to the lack of long-term records for the purposes of validation, a similar system was designed to synthesized TWSC (i.e., TWCC since there is no groundwater and surface water storages) from the Noah driven GLDAS land surface model, which is endowed with long-term series in its Version 2. This network was trained to predict TWCC from 2000 to 2010, after which, it was used to backcast TWCC from 2010 to 1979. Records from 1979 to 1999 from the reconstructed and the original datasets were then used to validate results from the ANN NARX.

Overall, the network employed to reconstruct Noah series yielded good spatial accuracies, with the spatially-averaged TWCC estimates presenting median RMSE, NSE and $R^2$ of 8.06 mm/month, 0.76 and 0.88 respectively. Thus, the artificial neural network method proved robust to adequately reconstruct GRACE-derived TWS estimates over West Africa. For the real GRACE data, the reconstructed TWSC series, showed strong agreement with the original GRACE data, presenting similar trends over their overlapping periods. The results, upon validation showed high skill scores, mostly in areas south of the Sahara Desert, but low performance in the desert dues to the very low amplitudes in signals. The spatially-averaged series for major basins and sub-climatic zones, as well as, the whole of West Africa itself, presented median RMSE, NSE, and $R^2$ of 11.83 mm/month, 0.76 and 0.89, respectively. These results agree with those of the “closed-loop” simulation and thereby we can conclude that the NARX network method used here proved robust to adequately reconstruct GRACE-derived TWSC estimates over West Africa.