Quantification of surface water extent and volume in the Inner Niger Delta (IND) over 2000–2022 using multispectral imagery and radar altimetry

Abstract

Spatio-temporal dynamics of surface water reservoirs at regional and global scales remain poorly understood. Using satellite remote sensing provides a unique opportunity to address this problem. This study aims to (1) quantify the extent and volume of surface water and (2) compare our approach with other datasets. We utilized MODIS-based spectral indices to monitor temporal variations in inundation extent. By interpolating water levels across the surface water extent, we generated water level maps and quantified surface water volume from 2000 to 2022. Evaluation against ICESat-2 data involved 64 comparisons, with approximately 58% showing an R² value greater than or equal to 0.5, and 38% were higher than or equal to 0.7. Compared to the flooding model, our method aligns more closely with ICESat-2 data, contrary to the flooding model which tends to overestimate water levels and, consequently, water storage.

Keywords:

• Surface water extent
• surface water volume
• Inner Niger Delta
• multispectral imagery
• radar altimetry

1. Introduction

The surface water reservoir, consisting of lakes, rivers, floodplains, and wetlands, plays a significant role as one of the primary water resources for ecosystems and populations (Kundzewicz et al. 2007; Younger 2012; Calvin et al. 2023). Within these reservoirs, wetlands, including floodplains, have a major role in regulating river flows by storing large amounts of water for extended periods, sometimes lasting several months (Bourrel et al. 2009; Fassoni‐Andrade et al. 2021; Kreibich et al. 2022). They also contribute to air temperature modulation and increased evapotranspiration (Krinner 2003). From a geochemical perspective, spatial and temporal variations in floodplains are of great importance for methane emissions (Bousquet et al. 2006) and carbon trapping and release (Abril et al. 2014). However, the spatiotemporal dynamics of floodplains are still poorly understood at regional and global scales. The main reasons for this lack of understanding are the absence of multi-year measurements of flooded area and water height, as well as limited incorporation of these variables in hydrological and hydrodynamic models.

Satellite remote sensing offers the possibility of quantifying surface water stocks using various techniques (Papa and Frappart 2021; Cretaux et al. 2023). Besides the increase in frequency of extreme droughts and floods (e.g. Kreibich et al. 2022) exacerbates the importance of remotely sensed observations to spatialize the signature of such events. The use of SAR interferometry (Alsdorf et al. 2007; Lu and Kwoun 2009; Hong and Wdowinski 2011; Wdowinski and Hong 2015) allows for high-resolution mapping of height variations (<100 m) between two satellite passes but is limited by low temporal repeatability and small coverage areas. The other two techniques rely on combining flooded extents provided by satellite imagery with either digital elevation models (DEMs) or radar altimetry (e.g. (Kitambo et al. 2023). These approaches have primarily been employed using passive microwave-derived flooded extents GIEMS (Global Inundation Extent from Multi-Satellite) with a spatial resolution of 25 km (Prigent et al. 2020). Their main limitations are the low spatial resolution of flood maps and the limited temporal depth due to the use of altimetric data from a single mission most of the time (Frappart et al. 2012; Papa et al. 2015). Consequently, these different methods do not provide long time series of floodplain water storage dynamics. It is necessary to combine missions with suitable spatial and temporal resolutions to properly monitor the hydrological cycle of floodplains. Recent studies have created longer time series (since 2000) using MODIS multispectral images (spatial resolution of 500 m and temporal resolution of 8 days) and records from several altimetric missions such as ERS-2/ENVISAT/Saral. These studies have been conducted in the Arctic environment of the Mackenzie Delta between 2000 and 2015 (Normandin, Frappart, Lubac, et al. 2018), the Tonle Sap, largest lake of East Asia, between 1993 and 2016 (Frappart et al. 2018), the Mekong Basin between 2003 and 2010 (Pham-Duc et al. 2017), and more recently, between 2000 and 2020. However, these initial studies need to be generalized to other types of environments at different scales and over longer time series. Nevertheless, so far, this method has not been validated due to a lack of in situ or ancillary data. Therefore, the main objectives of this study are to (1) quantify surface water extent and volume and (2) compare/validate them. For this purpose, the study was carried out on the Inner Niger Delta, located in Mali, since a digital flooding model based on in situ data was available.

2. Materials and methods

The different steps of the method used in this study are presented in Figure 1. Firstly, surface water extent has been calculated using multispectral imagery and water levels using radar altimetry. Then, the combination of these two data have been interpolated to calculate surface water volumes time series over 2000–2022. Finally, all these data have been compared using existing methods (Figure 1).

Figure 1. Flowchart of the method used.

2.1. Study site

The Inner Niger Delta (IND), located in Central Mali (between 3–5° W and 13–17° N), is the second-largest wetland in Africa, and supporting over one million herders, fishermen and farmers (De Noray 2003). The IND is composed of a variety of swamps, rivers, numerous lakes and non-permanent flooded areas which formed an important hydrological network, localized in a flat and sandy basin at 250 meters of altitude (Figure 2). The IND is fed by the Niger River, entering by Ke-Macina in the delta (Figure 2), and its main tributary the Bani, entering by Sofara (Figure 2), all both covering a total drainage area of 73,000 km² (Thompson et al. 2021). Both Niger and Bani Rivers originate from the South, taking their sources on the high plateau of Guinea at the foot of the Tingi Mountains.

Figure 2. Inner Niger Delta (IND) study site and location of main cities. In situ station (black circles), ERS-2/Envisat/Saral missions (purple triangles), Sentinel-3A (green crosses) and Sentinel-3B (black crosses) are represented on the map. The Niger and Bani rivers are in blue. Flow direction is indicated by blue arrows. The bottom of the map is the digital elevation model Shuttle Radar Topography Mission (SRTM, freely available on NASA website).

The IND is under the affluence the West African Monsoon, at the origin of an excess of rainfall between June and September in the Bani and Niger rivers in the upstream, with annual total in average of 250 mm over this period (Thompson et al. 2021). Precipitation is seasonal and characterized by large inter-annual variability (Zwarts and Grigoras 2005; Thompson et al. 2021). This rainfall excess produces an average annual discharge (cumulated) at the delta’s entries of 1490 m³.s⁻¹ over 1955–1996 period (Bergé-Nguyen and Crétaux 2015). This is at the origin of a wet season with flood events from August to December with a maximum flood extent that can reach 30,000 km² (Zwarts and Grigoras 2005; Zwarts et al. 2006; Bergé-Nguyen and Crétaux 2015). The maximum flood extent strongly varies over time, depending on the amount of rainfall in the upstream river’s watersheds. During flood events, water levels strongly increase of 4.5–7 m compared to low waters. During the dry period from March to May, the flooded area dries out, except the river mainstem and permanent lakes. The flow is reduced to nearly nothing, on average 80 times less important than during rainy season (Davids et al. 2019).

2.2. Data collection and methods

All the satellite missions used in this study are presented in Table 1 and are detailed in the following subsections.

Table 1. Satellite products used in the study.

Mission	Terra	ERS-2 ENVISAT	Saral	Sentinel-3A and 3B	ICESat-2	Global Surface Water Dynamics (Pickens et al. 2020)	Global Surface Water (Pekel et al. 2016)
Instrument	MODIS	Radar Altimeter (RA) Radar Altimeter (RA-2)	AltiKa	Sar Radar Altimeter (SRAL)	Advanced Topographic Laser Altimeter System (ATLAS)	Landsat-5, 7 and 8 missions	Landsat missions
Product	MOD09A1 composite Level-3	GDR Level-2	GDR Level-2	GDR Level-2	ATLAS/ICESat-2 L3A ATL13	GSWD	GSW
Downloading	https://appeears. earthdatacloud. nasa.gov/	https://hydroweb. next.theia-land.fr/	https://hydroweb. next.theia-land.fr/	https://hydroweb. next.theia-land.fr/	https://nsidc.org/ data/icesat-2	https://glad.umd.edu/ dataset/global-surface- water-dynamics	https://global- surface-water.appspot. com/download
Data type	Multispectral imagery	Radar altimetry	Radar altimetry	Radar altimetry	Lidar	Multispectral imagery	Multispectral imagery
Use objective	Flood map	Water level	Water level	Water level	Validation	Validation	Validation
Operation date	2000–now	1995–2003 2002–2012	Since 2013	2016–now 2018–now	2018–now	1999–2021	1984–2021
Space agency	NASA	ESA	CNES, ISRO	ESA	NASA	–	–
Repetitivity (days)	8	35	35	27	91	30	30
Spatial reoslution (m)	500	–	–	–	–	30	30
Spectral indexes	EVI = $a\times \frac{\rho \mathit{\text{NIR}}–\rho \mathit{\text{red}}}{\rho \mathit{\text{NIR}}+b\times \rho \mathit{\text{red}}-c\times \rho \mathit{\text{blue}}+d}$ LSWI = $\frac{\rho \mathit{\text{NIR}}–\rho \mathit{\text{SWIR}}}{\rho \mathit{\text{NIR}}+\rho \mathit{\text{SWIR}}}$	–	–	–	–	NDWI = $\frac{\rho \mathit{\text{NIR}}-\rho \mathit{\text{green}}}{\rho \mathit{\text{NIR}}+\rho \mathit{\text{green}}}$ MNDWI = $\frac{\rho \mathit{\text{green}}-\rho \mathit{\text{SWIR}}}{\rho \mathit{\text{green}}+\rho \mathit{\text{SWIR}}}$	NDVI = $\frac{\rho \mathit{\text{red}}–\rho \mathit{\text{NIR}}}{\rho \mathit{\text{NIR}}+\rho \mathit{\text{SWIR}}}$ Hue-Saturation- Value (HSV)

2.2.1. Surface water extent from multispectral imagery

2.2.1.1. MODIS product

The Moderate Resolution Imaging Spectroradiometer (MODIS) sensor is part of the payload of Terra and Aqua satellites launched in 1999 and 2002 respectively. Here, we used the MOD09A1 version 6.1 product of level 3 (Vermote 2015), freely available at https://appeears.earthdatacloud.nasa.gov/. In the website AppEEARS, we uploaded our shapefile as input and a reprojection was made in the projection chosen to merge different tiles corresponding to the study site. The product obtained in NetCDF is level-3 product, corrected from atmospheric conditions such as gasses, aerosols and Rayleigh scattering. Surface spectral reflectance of 7 spectral bands, from visible (blue: 459–479 nm, green: 545–565 nm, red: 620–670) to Short Wave Infra-Red (SWIR1: 1628–1652 nm, SWIR2: 2105–2155 nm), are provided in the product. For each pixel of a 500-m spatial resolution, best surface reflectance is selected from all acquisitions within the 8-day composite-period, based on a low view angle, the absence of clouds or cloud shadow and aerosol loading. For this study, a total of 1028 composites have been used between 2000 and 2022 over the IND.

2.2.1.2. Flood map method

Different methods have been developed in recent years for the detection of water surfaces, mainly based on active (Radar and LIDAR) or passive (multispectral and hyperspectral imagery) sensors (Albertini et al. 2022; Sogno et al. 2022). Here, we used multispectral imagery to detect surface water. The method used is based on the one developed by Sakamoto et al. (2007) and adapted by Normandin, Frappart, Lubac et al. (2018) in the Mackenzie Delta, Tonle Sap Lake (Frappart et al. 2018) and more recently in the Mekong River basin. Thresholds on the Enhanced Vegetation Index (EVI), Land Surface Water Index (LSWI) and the difference value between EVI and LSWI permit to classify each pixel in a category (non-flooded, mixed, flooded and permanent water body). More details are shown in Figure 1 and in Sakamoto et al. (2007) and (Normandin, Frappart, Lubac et al. 2018). The two spectral indexes are defined as follow (Huete 1997; Xiao et al. 2005): (1) $$\mathrm{\text{LSWI}}=\frac{\rho \mathit{\text{NIR}}‐\rho \mathit{\text{SWIR}}}{\rho \mathit{\text{NIR}}+\rho \mathit{\text{SWIR}}}$$(1) (2) $$\mathrm{\text{EVI}}=a\times \frac{\rho \mathit{\text{NIR}}‐\rho \mathit{\text{red}}}{\rho \mathit{\text{NIR}}+b\times \rho \mathit{\text{red}}-c\times \rho \mathit{\text{blue}}+d}$$(2)

Where:

• $\rho \mathit{\text{NIR}}$ is the surface reflectance value in the NIR (841–875 nm, band 2)

• $\rho$SWIR is the surface reflectance in the SWIR (1628–1652 nm, band 6)

• $\rho$ red is the surface reflectance value in the red (621–670 nm, band 1)

• $\rho$ blue is the surface reflectance value in the blue (459–479 nm, band 3)

• a, b, c and d are equal to 2.5, 6, 7.5 and 1, respectively.

2.2.2. Water levels from radar altimetry

Originally developed to study ocean topography, radar altimetry has been rapidly used for studying continental waters for over two decades now (Birkett 1995; Frappart et al. 2006; Crétaux et al. 2011; Da Silva et al. 2012). The measurement of water levels is based on the measurement of the round-trip time of electromagnetic waves sent by the radar altimeter onboard the satellite. This measurement allows to estimate the altimeter range which corresponds to the distance between the satellite and the rivers, lakes, and floodplains in our case (Chelton et al. 2001; Escudier 2017). By applying atmospheric and geophysical corrections to the range, it is then possible to deduce the height of water levels as the difference between sattelite orbit and the corrected range. Detailed explanations of the principle can be found in Crétaux et al. (2017).

Radar altimetry has been used to obtain long time series of water level data since the early 1990s. Virtual Stations (VS) of water level are created by selecting the intersections between the altimeter tracks and the objects of interest. In this study, we utilized the 27-day and 35-day missions. The VS from the 35-day missions are from the ERS-2/Envisat/Saral missions spanning from 1995 to 2016, as created in the study by (Normandin, Frappart, Diepkilé et al. 2018). They amount to 52/63/62, respectively. The 27-day VS are derived from the Sentinel-3A and 3B missions, covering the period from 2016 to the present, and obtained from the study by Diepkilé et al. (2021). There are 28 VS for Sentinel-3A and 23 for Sentinel-3B.

The VS for the 27-day missions were generated using the AlTiS software (Frappart et al. 2021a), while the 35-day ones were created using the MAPS software (Frappart et al. 2015; Normandin, Frappart, Diepkilé et al. 2018). They were manually created in three main steps: a coarse delineation of the VS using Google Earth, visual inspection for a finer selection of valid elevation data, and computation of the water level time series (Normandin, Frappart, Diepkilé et al. 2018). Surface water volume combining multispectral imagery and radar altimetry

Figure 3. Time series of surface water extent (km²) calculated using NORMANDIN2018, ZWARTS 2018, PEKEL2016 and PICKENS2020.

Using surface water extent and water levels, maps of surface water levels were obtained. The method was developed by (Frappart et al. 2005, 2012) and consists in interpolating altimetry-based water stages over inundates using an inverse distance weighting technique. This method has already been applied in different study sites, as in the Negro (Frappart et al. 2005), and Amazon (Frappart et al. 2012) river basins, the Mackenzie delta (Normandin, Frappart, Lubac et al. 2018), Mekong basin. Finally, time series of surface water volumes (km³) was calculated over 2000–2022 in the IND using the following equation (Frappart et al. 2012; Normandin, Frappart, Lubac et al. 2018): (3) $$\text{Δ}V={\displaystyle \sum _{jϵs}\left[\text{h}\left({\lambda }_j,{\phi }_j\right)-h_{\text{min}}\left({\lambda }_j,{\phi }_j\right)\right].}\hspace{0.17em}{\delta }_j.\Delta S$$(3)

With:

• $\mathrm{\Delta }V:$ anomaly of surface water volume (km³)

• $S:$ is the surface of the Inner Niger Delta (km²)

• h(λ_j, φ_j): the water level at longitude λ and latitude φ

• h_min(λ_j, φ_j): the minimal water level for the pixel of coordinates (λ_j, φ_j)

• δ_j: equals 1 if the j^th pixel is associated with inundated and 0 if not

• ΔS: pixel size expressed as 0.25 km²

2.2.4. Comparison/validation datasets

2.2.4.1. Global Surface Water

The Global Surface Water (GSW) product, provided by the European Commission’s Joint Research Centre (JRC) under the framework of the Copernicus Program, was used to compare surface water extent using (Pekel et al. 2014, 2016) method with those obtained by our method. The product contains six different layers: occurrence, occurrence change intensity, seasonality, recurrence, transitions and extent. Here, we only focused on the monthly surface water extent maps, available from 1984 to 2021 (updated data), with a spatial resolution of 30 m. In total, 261 images were downloaded for free here: https://global-surface-water.appspot.com/download. The surface water extent detection method is based on the use of a wide range of spectra, the Normalized Difference Vegetation Index (NDVI), and Hue-Saturation-Value (HSV) color-space transformations for different band combinations (SWIR2, NIR, red, NIR/green/blue, (Pekel et al. 2014)).

2.2.4.2. Global Surface Water Dynamics 1999–2021

Global Surface Water Dynamics (GSWD) maps, provided by Global Land Analysis & Discovery (GLAD) laboratory, were used to compare with surface water extent from our method. Global maps derived from Landsat-5, 7 and 8 scenes were used between 1999 and 2021 to characterize inter-annual and intra-annual open surface water dynamics. Different classification models utilize all the Landsat bands and normalized difference ratios of each pairs of bands (Pickens et al. 2020). These ratios include water indices such as NDWI (McFEETERS 1996) and MNDWI (Xu 2006). A total of 264 monthly percent water presence maps have been used between 2000 and 2021 at a spatial resolution of 30 m. All data are freely available at: https://glad.umd.edu/dataset/global-surface-water-dynamics.

2.2.4.3. ICESat-2

ICESat-2 (Ice, Cloud and land Elevation Satellite 2), part of the Earth Observating System (EOS) of NASA, is a mission for measuring ice sheet thickness, topography, vegetation characteristics and clouds (Markus et al. 2017). This satellite was launched on September 2018 by NASA and has a 91-day repeat cycle. It carries aboard the Advanced Topographic Laser Altimeter System (ATLAS), which is a space-based lidar and measures the travel time of laser photons from the satellite to the Earth and back (Markus et al. 2017). ATLAS is a low energy multibeam laser operating at a wavelength of 532 nm (green) and generates six beams arranged in three pairs to provide more ground coverage identified by its orientation (Left-L or right-R) and its spot number (1–3): GT1L, GT1R, GT2L, GT2R, GT3L, GT3R. In this study, we used the ATLAS/ICESat-2 L3A ATL13 product, which contains along-track surface water products for inland water bodies (Jasinski et al. 2021). Inland water bodies include lakes, reservoirs, rivers, estuaries and a 7 km near-shore buffer. In this product, elevation was provided with reference to WGS84. Data is freely available at: https://nsidc.org/data/atl13/versions/2#anchor-1.

2.2.4.4. Flooding model

Zwarts et al. (Zwarts and Grigoras 2005) initially used topographic maps to generate a flood map, which was further refined by employing Landsat TM satellite images with a resolution of 30 meters to create hydrographic maps. The identification of land and water was accomplished by analyzing spectral bands 5 (Near Infrared (NIR), 0.85–0.88 µm) and 7 (Shortwave Infrared 1 (SWIR 1), 2.11–2.29 µm). Specifically, a pixel was classified as water if the reflectance values fell within the range of 0.1–0.135 for band 5 and 0.7–0.9 for band 7; otherwise, it was considered as land. These maps were then combined to establish the connection between flood extent and water level using in situ records from Akka gauge station from 1992 to 2018. The measured water heights were provided by the Direction Nationale de l’Hydraulique (DNH) of Mali through in situ monitoring. By integrating all of these maps, a flooding model was constructed (Zwarts and Grigoras 2005; Zwarts et al. 2006). Subsequently, a flood prediction tool was developed based on this model, and it can be accessed freely at https://www.opidin.org/fr. This tool serves various purposes, including functioning as a warning system for communities residing in the IND. More recent advancements have been made to the flooding model, including the expansion of the model coverage area, utilization of a larger dataset, and incorporation of additional in situ water level stations (Davids et al. 2019).

3. Results

3.1. Surface water extent

Water surface extents for the period 2000–2022 were computed using the approach initially developed by Sakamoto et al. (2007) and later adapted by Normandin, Frappart, Lubac et al. (2018), the method proposed by Pekel et al. (2014), and the method introduced by Pickens et al. (2020) (Figure 3). However, owing to the absence of in situ water level data measured at the Akka station between 2018 and 2022, the water surfaces calculated using the method by Zwarts and Grigoras (2005) and Davids et al. (2019) are limited to calculations up to 2018. To enhance clarity, the method developed by Sakamoto et al. (2007) and adapted by Normandin, Frappart, Lubac et al. (2018) is denoted as NORMANDIN2018, the method of Pekel et al. (2016) as PEKEL2016, the method of Pickens et al. (2020) as PICKENS2020, and the method by Zwarts et al. (2006) and Davids et al. (2019) as ZWARTS2018.

The water surface areas obtained by NORMANDIN2018, which considers flooded and mixed pixels, and by ZWARTS2018 exhibit a relatively similar pattern, with a mean difference of 19% over 2000–2018. Over the common period of availability of the two datasets (2000–2018), the maximum flooding peak is 17,026 km² in 2003 for NORMANDIN2018 and 18,634 km² in 2010 for ZWARTS2018, while the minimum flooding peak is 6,987 km² in 2018 for ZWARTS2018 and 5919 km² for NORMANDIN2018 (Figure 3). The extents obtained by PEKEL2016, NORMANDIN2018 considering only flooded pixels and PICKENS2020 are similar over the period 2000–2013, with maximum values in 2004 (7777 km², 4942 km² and 6457 km², respectively) and minimum in 2011 (∼3000 km² for the 3 methods, Figure 3). From 2013 onwards, the extents determined by PICKENS2020 are similar to those obtained by NORMANDIN2018 by considering flooded and mixed pixels.

The mean annual flood durations are shown in Figure 4. Figure 4(a) shows the mean annual flood duration obtained using the method of NORMANDIN2018 without considering mixed pixels, Figure 3(b) with mixed pixels, Figure 4(c) the method of ZWARTS2018, Figure 4(d) with the method of PEKEL2016 and Figure 4(e) the method of Pickens et al. (2020). In order to have an identical comparison period for all the methods, the average was taken over the years 2013–2018. All the methods highlight the main river, the Niger, and its main tributary, the Bani, shown in dark blue with differences in flood duration between upstream and downstream for NORMANDIN2018, shorter durations for PEKEL2016 and longer durations for PICKENS2020 on rivers. Lakes are also identified. The main differences are found at the level of the floodplain, located in the central part: although the method of ZWARTS2018 provides larger estimates of the extent of the flood compared to the other methods, and especially the method of NORMANDIN2018 without the mixed pixels which provides lower flood extent. Furthermore, NORMANDIN2018 does not identify the network as a whole, and in particular all the tributaries present in latitudes 13–15°N, which are visible using the methods of PEKEL2016 and PICKENS2020.

Figure 4. Mean flooding duration using the five methods. (a) from NORMANDIN2018 (with mixed pixels, (b) from NORMANDIN2018 (with mixed pixels), (c) from ZWARTS2018, (d) from PEKEL2016, and (e) from PICKENS2020.

3.2. Water level maps

3.2.1. Comparison between water level maps from NORMANDIN2018 and ICESat-2

The water level maps produced by NORMANDIN2018, obtained combining the surface water extents with the water level time series obtained from radar altimetry, were compared with those obtained from ICESat-2 data. Each time the ICESat-2 satellite passes, a transect of levelled water levels is obtained. By retrieving the coordinates from NORMANDIN2018 water level maps closest to those of ICESat-2, it is then possible to calculate statistics for each transect: the statistical parameters number of data (N), coefficient of determination (R²), root mean square error (RMSE) and bias (if the data are levelled) were calculated and are presented in the histograms in Figure 5(a–d), respectively. As a consequence, comparisons between our method and the flooding model, then the flooding model and ICESat-2 data, only the bias is not estimated as the flooding model is not levelled to an ellipsoid or a geoid model. We note that the comparison between our method and the ICESat-2 data presents the results with the highest R². Indeed, of the 64 comparisons, 29 have an R² determination coefficient greater than 0.6 and 18 between 0.2 and 0.6. In addition, 45 comparisons had a bias of between −0.5 and 0.5 m. Finally, 38 comparisons have a RMSE value between 0.25 and 0.75 m.

Figure 5. Comparison of water level maps from NORMANDIN2018 and ICESat-2, ZWARTS2018 and ICESat-2. For each comparison, (a) N is the number of samples compared averaged used for each comparison, (b) is the determination coefficient R², (c) is the Root mean Square Error (RMSE) in meters and (d) is the bias in meters.

Comparison of water levels measured by ICESat-2 and NORMANDIN2018 along altimetry tracks are shown for two cases (Figure 6), one in the downstream part (Figure 6(b)) and another in the upstream part (Figure 6(c)). For each example, water levels from NORMANDIN2018 show similarities with the ones obtained using ICESat-2, especially for the comparison in the upstream part, with a R² = 0.978 and a bias of 0.138 m.

Figure 6. Comparison of water levels measured (m) by ICESat-2, NORMANDIN2018 and ZWARTS2018 along two altimeter tracks, (a) map of the locations of the two altimetry tracks chosen, (b) water levels according to the latitude in the downstream and (c) according to the upstream. For each figure, ICESat-2 points are shown in blue, NORMANDIN2018 in green and ZWARTS2018 in black. For each comparison between ICESat-2 and NORMANDIN2018, N, R² and bias are displayed in green and in black for ICESat-2 and ZWARTS2018 comparison.

The spatialization of the statistical parameters between our method and ICESat-2 are presented in Figure 7. Figure 7(a) corresponds to the number of comparisons N, Figure 6(b) is the R², Figure 7(c) is the RMSE (in meters) and Figure 7(d) is the bias (in meters). The coefficient of determination R² is lowest in the downstream section, just outside the delta (Figure 7(b)). The highest bias values can be seen in the upstream part of the zone (Figure 7(d)). The same applies to the RMSE (Figure 7(c)).

Figure 7. Maps of statistical parameters for the water level maps from NORMANDIN2018 (combining surface water extent and radar altimetry) with ICESat-2 data, superimposed on a Google Earth image from ESRI world imagery in QGIS software. (a) N is the number of samples compared, (b) the determination coefficient R², (c) the Root mean Square Error (RMSE, in meters) and (d) the bias (in meters).

3.2.2. Comparison between water level maps from ZWARTS2018 and ICESat-2

Comparisons were made between ICESat-2 water levels and ZWARTS2018 inundation model. Due to the lack of in situ data after 2018, only 10 comparisons between in situ data and models could be achieved. The results are presented in Figure 4(a,b) for the number of samples N for each comparison and the associated determination coefficient R², respectively. Despite the high number of samples obtained for each comparison (with a median value of 258), the determination coefficients obtained are less than 0.1.

3.3. Surface water volume

Time series of surface water volumes were calculated using ZWARTS2018, NORMANDIN2018, PEKEL2016 and PICKENS2020 and are presented in Figure 8. As in situ water level data measured at Akka are not available after 2018, ZWARTS2018 time series stops at this period as these data are required to run the inundation model. The dynamics in surface water volumes are higher using ZWARTS2018 than our method, with an average difference of +21%. PICKENS2020 surface water volume are quite similar at those obtained using NORMANDIN2018 since 2013, and before they are similar to those of PEKEL2016. In all methods, intra- and inter-annual variations are observed, with maximum values found in November. The variations are the same as those already observed in the time series of water surfaces, shown in Figure 3. Over the period 2000-2022, the maximum volume observed in 2018 is equal to 34 km³ according to ZWARTS2018, 26 km³ using NORMANDIN2018, 22 km³ using PEKEL2016 and 17 km³ using PICKENS2020. The minimum volume in 2011 is estimated at 14 km³ by ZWARTS2018, 8 km³ by NORMANDIN2018, 0.5 km³ using PEKEL2016 and 0.3 km³ using PICKENS2020 (Figure 8).

Figure 8. Time series of surface water volume (km³) calculated using NORMANDIN2018, ZWARTS2018 PEKEL2016 and PICKENS2020 methods.

4. Discussion

Surface water extents were quantified using the method of NORMANDIN2018 (taking wet pixels on one side, and wet and mixed pixels on the other), the flooding model of ZWARTS2018, PEKEL2016 and PICKENS2020. The average annual flood duration over the period 2013–2018 was calculated using the different methods, and some differences stand out. Firstly, using our method, considering the mixed pixels makes it possible to better identify the floodplains located in the centre of the Inner Niger Delta (IND), between latitudes 14–15°N and longitudes 4–5°W. The patterns obtained are relatively similar to those obtained by PEKEL2016 and PICKENS2020, but our method differs from the other ones in the upstream of the IND, since the main and secondary rivers have a lower flood duration than the rest of the river. This observation can be explained by the fact that at this downstream section, the width of the Niger River is less than 500 m, whereas in the center of the delta there are large floodplains of several kilometers wide. Furthermore, due to the 500 m spatial resolution of the MODIS data used in NORMANDIN2018, the secondary network is not as well identified as by the methods of PEKEL2016 and PICKENS2020, which have both a spatial resolution of 30 m. PEKEL2016 product shows lower flood duration in the main river, compared to PICKENS2020 and ZWARTS2018 products. These differences are likely to be due to the method used to map the extent of water surfaces. PICKENS2020 method uses the NDWI spectral index, and ZWARTS2018 uses spectral bands (NIR and SWIR) related to water, and taking less into account the presence of water under the vegetation. The methods of NORMANDIN2018 and PEKEL2016 considers other spectral indices that also accounts for the presence of vegetation, making possible to better discriminate between water and vegetation pixels. Thus, the methods of PICKENS2020 and ZWARTS2018 provides larger estimates of water surfaces in the IND than NORMANDIN2018 and PEKEL2016. Before 2013, surface water extent from PICKENS2020 are lower due to the presence of stripes on Landsat-7 images, reducing the extent of water surfaces. Much larger estimates are obtained with the flooding model of ZWARTS2018 than for the other methods.

Subsequently, water level maps produced by NORMANDIN2018 and the one from ZWARTS2018, were compared to ICESat-2 along-track height profiles (Figure 5). The results of the comparison between our method and ICESat-2 show a histogram with higher R² values, 45% with higher values than 0.6 (Figure 5(b), in blue). Out of the three data sources, the comparisons between ZWARTS2018 and ICESat-2 shows the lowest R² values (Figure 5(b), in green), with an R² less than 0.2 for all the cross-sections. Water levels along ICESAT-2 altimetry tracks (Figure 6) show larger differences between ICESat-2 and ZWARTS2018 compared to ICESat-2 and NORMANDIN2018 in the downstream and upstream parts of the IND. The discrepancies observed ZWARTS2018 and ICESat-2 along-track water profile are most likely due to the use of a unique water level records to derive the water elevation and storage in the whole IND. Due to the propagation time of the flood pulse, estimated between 20 and 80 days depending on the year (Mariko et al. 2013), there is a time-lag in the different phases of the flood pulse (rise, maximum, decrease) between the upstream and the downstream parts of the IND.

Changes in surface water volumes were calculated using ZWARTS2018 over the period 2000–2018 and NORMANDIN2018 over the period 2000–2022. PEKEL2016 and PICKENS2020 were used to quantify changes in surface water volumes from 2000 to 2021. The changes in surface water volumes obtained using ZWARTS2018 method are 21% higher over the common period 2000–2018 than those calculated using NORMANDIN2018. The difference could be attributed to the absence of a digital elevation model (DEM) in the mapping of water pixels, leading to the exclusion of slope consideration in ZWARTS2018. Despite a relatively flat area, the slope between the entrance and exit of the IND is estimated at 0.2%. As previously discussed, this higher estimation of surface water volumes is most likely related to the hydrodynamics of the delta (i.e. the propagation of the flood in the IND). Changes in surface water volumes from PICKENS2020 and PEKEL2016 are lower to those obtained using NORMANDIN2018 and ZWARTS2018, approximately divided by two. These differences are probably due to the fact that Landsat-7 images were used in these methods and present stripes over images, strongly reducing the surface water extent before 2013. From 2000 to 2013, the average difference between the surface water volumes calculated by NORMANDIN2018 and those calculated by PEKEL2016 is +147%, and +137% compared to PICKENS2020. Over the period 2014–2021, the average difference is +5% between NORMANDIN2018 and PEKEL2016, and +24% between NORMANDIN2018 and PICKENS2020. Thus, despite the fact that our NORMANDIN2018 method is at a spatial resolution of 500 m, the differences with the Landsat data, at a spatial resolution of 30 m, from PEKEL2016 and PICKENS2020 are small from 2014 onwards. In addition, our method provides data since 2000, unlike the Landsat data, which has a problem before 2013.

The dynamics of filling and emptying of the IND were observed through variations in the water surface extent and changes in surface water volumes. Mean seasonal cycle between surface water extent (km²) and surface water volume (km³) using NORMANDIN2018 and ZWARTS2018 over 2000–2018 are displayed in Figure 9 showing two different hydrological regimes associated to the two methods, with a cycle for NORMANDIN2018 (Figure 9(a)) and rather a linear pattern for ZWARTS2018 (Figure 9(b)). The ZWARTS2018 method pattern shows some slope discontinuities compared to NORMANDIN2018. Regarding NORMANDIN2018 (Figure 9(a)), the mean annual cycle can be divided into four phases: a lower phase from February to May, a rising phase from June to August, a high-water phase in September and a rising phase from October to January. Compared to ZWARTS2018, the main difference is observed for the maximum phase, which occurs in October instead of September for NORMANDIN2018. Thus, the water filling duration in the Inner Niger Delta seems to be faster for NORMANDIN2018 (3 months, Figure 9(a)) than for ZWARTS2018 (4 months, Figure 9(b)) and the water draining duration is faster in ZWARTS2018 (3 months, Figure 9(b)) than in NORMANDIN2018 (4 months, Figure 9(a)). Accordingly, the hydrological regimes will be different depending on the method used to quantify the surface water extent and volume of water.

Figure 9. Mean seasonal cycle between surface water extent (km²) and surface water volume (km³) using (a) NORMANDIN2018 and (b) ZWARTS2018. The numbers represent the month of the year.

NORMANDIN2018 method shows some advantages to use compared to other. Firstly, the method makes it possible to obtain long time series over approximately 25 years, and thus to carry out studies on the impacts of climate change and human activities. No in situ data is required, contrary to the hydrological models. In addition, the NORMANDIN2018 method is relatively simple to implement, since the surface water mapping method is based on the use of spectral indices and thresholds. For the creation of temporal series of water levels, tools have been developed that make the processing of these data accessible. Finally, all the data used is accessible to everyone, free of charge

5. Conclusion

The temporal series of surface areas and surface water volumes were quantified using MODIS multispectral imagery (spatial resolution of 500 m, temporal sampling of 8 days) and radar altimetry (temporal sampling of 27–35 days depending on the mission) over the period 2000–2022. Initially, the water surfaces were compared with various existing products (PEKEL2016, ZWARTS2018 and PICKENS2020). The different methods show similar seasonal cycles (with maximum flooding in October/November) but different spatial patterns. Indeed, our method seems to have differences in flood duration between upstream and downstream of the Inner Niger Delta, with shorter flood durations observed by PEKEL2016 and longer by PICKENS2020. However, all of these methods make it possible to identify the main network of the Niger River and the tributaries, including the main one, the Bani. Next, the water level maps produced by our method (multispectral imagery and radar altimetry) were compared to the ICESat-2 data. Similar comparisons to the ICESat-2 along-track height profiles were achieved with ZWARTS2018. Comparison between NORMANDIN2018 and ICESat-2 had the highest R², with 45% of comparisons superior to 0.6 and a RMSE value between 0.25 and 0.75 m. The comparisons between ZWARTS2018 and ICESat-2 data show the lowest R², below 0.2.

Therefore, for the first time, NORMANDIN2018 has been compared to different types of data: PEKEL2016, PICKENS2020, and ZWARTS2018. NORMANDIN2018 combining multispectral imagery and radar altimetry has been validated. Following this study, this method can be applied to other large basins or generalised. Water surfaces can be calculated using the mapping method, if the study site does not have to much clouds and vegetation, then combined with altimetric water heights to obtain time series of surface water volumes. All these data (MODIS multi-spectral imagery and radar altimetry) are available free of charge.

However, some improvements can be made at different levels:

• Surface water extent from multispectral imagery:

  • Improved spatial resolution by reducing the pixel size from 500 m to 250 m, and improved temporal resolution by combining MODIS data from the Terra satellite with data from the MODIS sensor on the Aqua satellite. In this way, twice as much data will be obtained, making it possible to obtain more data in areas where the presence of clouds is very high.

  • Data continuity: application of the surface water mapping method to the Sentinel-3 mission with the OLCI (Ocean and Land Colour Instrument) sensor at a spatial resolution of 300 m, with the aim of making the data sustainable in order to continue obtaining long time series when the Terra and Aqua satellites (MODIS sensor) will no longer be operational.

  • Smaller-scale study area: it would be possible to apply the method to higher spatial resolution images (such as Landsat-8/9 and Sentinel-2) to be able to work on smaller areas and at a finer scale (10/20/20 m).

• Water levels from radar altimetry:

  • Automatically generation of the creation of dense networks of virtual altimetric water level stations. The method was initiated by Frappart et al. 2021b and Enguehard et al. 2023, and would then make it easier to globalise the NORMANDIN2018 method.

• Surface water volumes:

  • Comparison with data from the SWOT (Surface Water and Ocean Topography) mission, launched in December 2022, with the aim of mapping continental surface water and estimating variations in surface water storage.

Disclosure statement

No potential conflict of interest was reported by the author(s).