Small Reservoirs, Landscape Changes and Water Quality in Sub-Saharan West Africa

: Small reservoirs (SRs) are essential water storage infrastructures for rural populations of Sub-Saharan West Africa. In recent years, rapid population increase has resulted in unprecedented land use and land cover (LULC) changes. Our study documents the impacts of such changes on the water quality of SRs in Burkina Faso. Multi-temporal Landsat images were analyzed to determine LULC evolutions at various scales between 2002 and 2014. Population densities were calculated from downloaded 2014 population data. In situ water samples collected in 2004 / 5 and 2014 from selected SRs were analyzed for Suspended Particulate Matter (SPM) loads, an integrative proxy for water quality. The expansion of crop and artiﬁcial areas at the expense of natural covers controlled LULC changes over the period. We found a very signiﬁcant correlation between SPM loads and population densities calculated at a watershed scale. A general increase between the two sampling dates in the inorganic component of SPM loads, concomitant with a clear expansion of cropland areas at a local scale, was evidenced. Results of the study suggest that two complementary but independent indicators (i.e., LULC changes within 5-km bu ﬀ er areas around SRs and demographic changes at watershed scale), relevantly reﬂected the nature and intensity of overall pressures exerted by humans on their environment, and locally on aquatic ecosystems. Recommendations related to the re-greening of peripheral areas around SRs in order to protect water bodies are suggested. The results reveal a strong negative correlation between cropland and sparse vegetation, and between dense vegetation and artiﬁcial area at the local scale (5-km bu ﬀ er areas). Dense vegetation and artiﬁcial surfaces were negatively correlated at the watershed scale as well, with also a negative correlation between dense vegetation and water at the same scale. Some expected correlations were not realized, e.g., negative correlation between cropland and sparse vegetation areas ( p = 0.023) at the watershed scale, and between sparse and dense vegetation areas ( p > 0.1) at both scales. The reduction in dense vegetation areas appeared only to be associated with an increase in artiﬁcial areas at both local and watershed scales. The confusion between cropland and artiﬁcial surfaces at the classiﬁcation


Introduction
Small reservoirs (SRs) are important artificial water resource structures in semi-arid environments around the world and noticeably in Sub-Saharan Africa [1]. They harvest rainwater during the rainy season and are often the only source of water for multiple uses such as domestic water supply, livestock watering, irrigation, etc., and also in delivering local employment opportunities during the dry season. In Burkina Faso, there are an estimated 1450 SRs (ca. 1/3 perennial), most of which were constructed between 1974 and 1987 in response to the Sahel droughts of the 1970s and 1980s [2]. These SRs mitigate adverse effects of rainfall variability [3], but may sometimes cause deleterious health effects [4][5][6]. 1.
What have been the dominant LULC changes occurring around SRs in Burkina Faso between 2002 and 2014? 2.
What are the possible linkages between LULC and demographic situations, and the water quality assessed in 2014 within a set of selected SRs? 3.
What have been the impacts of LULC changes observed between 2002 and 2014 on the water quality of a sub-series of these SRs studied twice, in 2004/5 and in 2014?
The first part of the paper deciphers LULC changes. The second one analyzes SPM measurements and correlates them with anthropogenic pressures and changes. Results are then discussed and recommendations are finally provided.

Materials and Methods
Study area: Burkina Faso (≈274,000 km 2 ) is a Sub-Saharan country located in West Africa ( Figure 1). Its population was estimated at 18.6 million in 2016, with an annual growth rate of 3.1% [54]. Remote Sensing Data: A total of twenty-two Landsat images (Enhanced Thematic Mapper, ETM+, and Operational Land Imager, OLI) from five scenes were downloaded from the Global Visualization Viewer (GloVis) and analyzed to map LULC changes between 2002 and 2014 (eleven images in both cases; Figure 2). In processing the data, only the reflective spectral bands (i.e., red, blue, green, NIR, SWIR1 and SWIR2) were extracted and analyzed. The raw digital numbers (DNs) were converted to at-sensor reflectance using the information in the associated metadata files to reduce the effects of atmospheric attenuation [59]. Reference data for classifying the Landsat images were primarily generated by interpreting various image color composites. The use of color composites (i.e., displaying different spectral bands as red, blue and green) permits the identification of general LULC classes such as the ones considered in this study with reasonable accuracy [60]. In instances where it was difficult to decipher one LULC type from the other (e.g., sparse from dense vegetation), Google Earth images were used. Five classes were considered for this study: (1) cropland, i.e., all cultivated land, (2) sparse vegetation, i.e., short grass with significant base soil areas, (3) dense vegetation, including woodland, grass, shrub and their mixture, (4) artificial areas, i.e., bares areas, settlements and (5) water bodies. Varying numbers of polygons ranging from 164 to 253 were generated for the Remote Sensing Data: A total of twenty-two Landsat images (Enhanced Thematic Mapper, ETM+, and Operational Land Imager, OLI) from five scenes were downloaded from the Global Visualization Viewer (GloVis) and analyzed to map LULC changes between 2002 and 2014 (eleven images in both cases; Figure 2). In processing the data, only the reflective spectral bands (i.e., red, blue, green, NIR, SWIR1 and SWIR2) were extracted and analyzed. The raw digital numbers (DNs) were converted to at-sensor reflectance using the information in the associated metadata files to reduce the effects of atmospheric attenuation [59]. Reference data for classifying the Landsat images were primarily generated by interpreting various image color composites. The use of color composites (i.e., displaying different spectral bands as red, blue and green) permits the identification of general LULC classes such as the ones considered in this study with reasonable accuracy [60]. In instances where it was difficult to decipher one LULC type from the other (e.g., sparse from dense vegetation), Google Earth images were used. Five classes were considered for this study: (1) cropland, i.e., all cultivated land, (2) sparse vegetation, i.e., short grass with significant base soil areas, (3) dense vegetation, including woodland, grass, shrub and their mixture, (4) artificial areas, i.e., bares areas, settlements and (5) water Water 2020, 12, 1967 4 of 25 bodies. Varying numbers of polygons ranging from 164 to 253 were generated for the respective classes, with an average of 189. Image classification was performed to distinguish the five general LULC classes. LULC classes have been assessed at two different geographical scales: within 5-km radius buffer areas around SRs, and at the watershed level (i.e., the upland-contributing area of each SR). 5-km radius buffer areas around SRs, and at the watershed level (i.e., the upland-contributing area of each SR).
Demographic data: Gridded population data of Burkina Faso were downloaded for the year 2014 from the WorldPop website [61]. The spatial resolution is of 100 m, and population data have been adjusted to United Nations population projection estimates. The accuracy of the data was assessed by comparing it with other widely used population datasets (e.g., GRUMP, GPW, LandScan, UNEP) based on data from four countries where official population estimates are reported at a higher administrative unit level than used in constructing the WorldPop dataset. The assessments revealed lower root mean square errors compared to the other datasets in all four countries [61]. Updated population estimates were extracted based on the 2014 gridded data. Population densities have been calculated at the two same geographical scales as those used for LULC classes' assessment. Water quality sampling: In 2014, 37 SRs, selected to represent the full range of sizes of water masses in Burkina Faso, were studied during the dry season (see [49] for a complete description). All the sampling sites are situated in the upper basin of the Volta River. Of 37 of these, 21 had already been studied in 2004/5 during the same season and in using the same sampling and analysis protocols ( Figure 1). Selected SRs were located close to main roads in order to minimize travel time between sample collection (one to three sites every day, early in the morning) and lab analyses (performed in the afternoon in Ouagadougou). One sample per site has been boat-collected during each campaign. Sampling locations were georeferenced with a Garmin GPS receiver ( Table 1). The depth of the water column at sampling locations was measured with a hand-held 0.1 m-resolution echo sounder (Plastimo Echotest II). Transparency was evaluated using a white Secchi disc (30-cm diameter), gently immersed and which depth of disappearance was measured (±2 cm). Temperature and conductivity vertical profiles were performed with a WTW LF 197 probe (accuracy ±0.5%) daily calibrated with an 84 µS/cm standard calibration solution. pH was also measured on the field with a WTW pH 197 (accuracy ±0.01 digit) daily checked against 4.01, 6.87 and 9.18 (Schott Instruments GmbH) buffer solutions. An AlgaeTorch ® (BBE-Moldaenke, GmbH) was used to quantify subsurface phytoplankton Demographic data: Gridded population data of Burkina Faso were downloaded for the year 2014 from the WorldPop website [61]. The spatial resolution is of 100 m, and population data have been adjusted to United Nations population projection estimates. The accuracy of the data was assessed by comparing it with other widely used population datasets (e.g., GRUMP, GPW, LandScan, UNEP) based on data from four countries where official population estimates are reported at a higher administrative unit level than used in constructing the WorldPop dataset. The assessments revealed lower root mean square errors compared to the other datasets in all four countries [61]. Updated population estimates were extracted based on the 2014 gridded data. Population densities have been calculated at the two same geographical scales as those used for LULC classes' assessment.
Water quality sampling: In 2014, 37 SRs, selected to represent the full range of sizes of water masses in Burkina Faso, were studied during the dry season (see [49] for a complete description). All the sampling sites are situated in the upper basin of the Volta River. Of 37 of these, 21 had already been studied in 2004/5 during the same season and in using the same sampling and analysis protocols ( Figure 1). Selected SRs were located close to main roads in order to minimize travel time between sample collection (one to three sites every day, early in the morning) and lab analyses (performed in the afternoon in Ouagadougou). One sample per site has been boat-collected during each campaign. Sampling locations were georeferenced with a Garmin GPS receiver ( Table 1). The depth of the water column at sampling locations was measured with a hand-held 0.1 m-resolution echo sounder (Plastimo Echotest II). Transparency was evaluated using a white Secchi disc (30-cm diameter), gently immersed and which depth of disappearance was measured (±2 cm). Temperature and conductivity vertical profiles were performed with a WTW LF 197 probe (accuracy ±0.5%) daily calibrated with an 84 µS/cm standard calibration solution. pH was also measured on the field with a WTW pH 197 (accuracy ±0.01 digit) daily checked against 4.01, 6.87 and 9.18 (Schott Instruments GmbH, Waltham, MA, USA) buffer solutions. An AlgaeTorch ® (BBE-Moldaenke, GmbH, Schwentinental, Germany) was used to quantify subsurface phytoplankton biomass, using chlorophyll a (Chl a, µg·L −1 ) as a proxy (see [49]). Then a 5-L subsurface grab water sample was collected for SPM concentrations measurement. The sample was filtered through a 200-µm mesh size net in order to eliminate large debris, then stored in a fridge (4 • C) and transported to the laboratory where it was processed within <1 to 5 h after collection. The SPM sampling approach described here applies to both 2004/5 and 2014 field campaigns.
LULC analysis: LULC dynamics between 2002 and 2014 around SRs and on their watersheds were analyzed by classifying Landsat data using the Random Forest (RF) classification algorithm [62]. RF was selected due to its superior performance over other classifiers in classification and regression problems [63,64]. The generated reference data/polygons were overlaid on the Landsat data and the corresponding pixel values (samples) extracted. For each LULC class, a maximum of 2000 samples were extracted. These samples were split into 50% training and 50% validation for the RF classification. The RF package [65] in the R statistical software [66] was used for the classification. RF generates multiple decision trees based on the training data, each tree built from a different sub-sample randomly selected from the training data. By creating multiple trees, RF reduces the correlation between trees in the "forest" and increases accuracy. In this study, five hundred decision trees were generated in each classification while other RF parameters were set to default values [62]. Although RF has an internal accuracy measure that is calculated from samples deliberately left out during tree building (i.e., the so-called 'out-of-bag' samples), the accuracy of each classification in this study was independently tested using the 50% validation samples. Overall, producer's and user's accuracies [67] were subsequently determined for each Landsat scene classified. LULC analysis was performed at two scales: watersheds and 5-km buffer areas around SRs.
Water quality analysis: Samples were collected to measure Suspended Particulate Matter (SPM) concentrations within SRs and to quantify their inorganic fraction (Particulate Inorganic Matter (PIM)). Back at the lab, two sub-samples of each specimen were filtered using Whatman GF/F filters (0.7 µm pore size), pre-weighed (0.1 mg) after being dried at 105 • C. Filtered volumes were adjusted in order to saturate filters without excess; they varied between 15 mL for Pouytenga (an urban site) and 700 mL for Bala (a natural reserve). After filtration, filters were dried again (105 • C for 90 min) and weighed to determine SPM concentration (mg·L −1 ). Filters were then placed in a furnace at 500 • C for four hours. They were weighed a second time for quantifying the loss on ignition and thus the inorganic fraction of SPM (PIM, expressed in percent of the SPM concentration). Duplicate measurements were averaged. The laboratory procedures described here (standard method ISO 11923:1997) applies to both 2004/5 and 2014 campaigns.
Statistical analysis: Spearman rank-order correlation coefficients were calculated to examine possible correlations between variables (LULC classes and population densities for the two geographical scales; environmental parameters and water quality proxies associated with the 37 reservoirs studied in 2014). Mann-Whitney non-parametric rank-sum tests were performed in order to compare changes in the water quality of SRs over time (2004/5 and 2014) and changes in anthropogenic pressures (LULC and demography) over the same period. The significance level was set at p < 0.01. Statistical analyses were performed using the software SigmaStat (v3.5, Systat Software Inc., San Jose, CA, USA).

LULC Changes (2002-2014) at Multiple Scales
Accuracy assessment of the LULC maps using the independent set of validation data revealed overall accuracies of more than 90% for six out of the 10 classifications and more than 81% for the remaining 4. Generally, classification accuracies were affected by confusion between sparse and dense vegetation on one hand, and cropland and artificial surfaces (especially bare areas) on the other. The latter case can be attributed to the use of images acquired during or after harvest (as was for this study). Indeed, cutting down the stalks of cereals (sorghum, maize) and transporting them for use elsewhere may leave the agricultural fields bare after harvest. This is aggravated as the dry season sets in. Confusion between cropland and vegetation types may also occur because medicinal, sacred, and fruit trees are maintained in croplands [15,68], which leads to mixed pixels. Table S1 (Supplementary materials) presents a summary of the classification accuracies obtained for all 10 classifications. Figure 3 shows box plots of LULC classes around the 37 SRs for 2002 and 2014. For both years, cropland and dense vegetation represent the dominant LULC classes. Cropland areas around SRs increased substantially between the two dates, while dense vegetation witnessed a significant reduction. Marginal but not statistically significant increases were observed for sparse vegetation, artificial surfaces and water bodies (p > 0.01). General trends seem globally persisting with very significant autocorrelations of all the different classes between the two dates (Spearman rank-order correlation coefficients varying between 0.478 and 0.829, p < 0.001). In other words, sites affected by important anthropogenic pressures in 2014 were already affected by important anthropogenic pressures in 2002.

LULC Changes (2002-2014) at Multiple Scales
Accuracy assessment of the LULC maps using the independent set of validation data revealed overall accuracies of more than 90% for six out of the 10 classifications and more than 81% for the remaining 4. Generally, classification accuracies were affected by confusion between sparse and dense vegetation on one hand, and cropland and artificial surfaces (especially bare areas) on the other. The latter case can be attributed to the use of images acquired during or after harvest (as was for this study). Indeed, cutting down the stalks of cereals (sorghum, maize) and transporting them for use elsewhere may leave the agricultural fields bare after harvest. This is aggravated as the dry season sets in. Confusion between cropland and vegetation types may also occur because medicinal, sacred, and fruit trees are maintained in croplands [15,68], which leads to mixed pixels. Table S1 (Supplementary materials) presents a summary of the classification accuracies obtained for all 10 classifications. Figure 3 shows box plots of LULC classes around the 37 SRs for 2002 and 2014. For both years, cropland and dense vegetation represent the dominant LULC classes. Cropland areas around SRs increased substantially between the two dates, while dense vegetation witnessed a significant reduction. Marginal but not statistically significant increases were observed for sparse vegetation, artificial surfaces and water bodies (p > 0.01). General trends seem globally persisting with very significant autocorrelations of all the different classes between the two dates (Spearman rankorder correlation coefficients varying between 0.478 and 0.829, p < 0.001). In other words, sites affected by important anthropogenic pressures in 2014 were already affected by important anthropogenic pressures in 2002.  LULC changes between 2002 and 2014 reported on an obvious trend toward the increase of anthropogenic pressures on landscapes around SRs. However, depending on the scale used for characterizing landscapes (i.e., 5-km buffer areas or watersheds), spatial autocorrelations between the different LULC classes differed (Table 2). Table 2. Spearman rank-order correlation coefficients between 2014 LULC classes at different scales around SRs (n = 37). Associated probabilities are indicated in italics and significant statistics in bold.

5-km Buffer
Sparse The results reveal a strong negative correlation between cropland and sparse vegetation, and between dense vegetation and artificial area at the local scale (5-km buffer areas). Dense vegetation and artificial surfaces were negatively correlated at the watershed scale as well, with also a negative correlation between dense vegetation and water at the same scale. Some expected correlations were not realized, e.g., negative correlation between cropland and sparse vegetation areas (p = 0.023) at the watershed scale, and between sparse and dense vegetation areas (p > 0.1) at both scales. The reduction in dense vegetation areas appeared only to be associated with an increase in artificial areas at both local and watershed scales. The confusion between cropland and artificial surfaces at the classification stage could have influenced this, as dense vegetation is most likely to be converted to cropland. However, increasing (historical) forest degradation and desertification, especially in the northern part of Burkina Faso [69], makes the negative correlation between dense vegetation and artificial surfaces putatively valid. The correlation between dense vegetation and water at watershed scale corroborates the existence of a link between the location of water masses and human pressures on watersheds (i.e., SRs seemed preferentially located in agriculture prone areas).

Indicators of Anthropogenic Pressures
Complementary and non-redundant indicators of population density and LULC changes that accurately represent anthropogenic pressures on SRs at local (5-km buffer area) and regional (watershed) scales were determined by analyzing correlations between the two datasets at both spatial scales. First, a comparison of 2014 population densities from the WorldPop data at watershed scales revealed a large diversity of demographic situations for the studied SRs (varying between <15 and >9000 Inh·km −2 - Table 1). Population densities calculated at local and watershed scales appeared highly correlated (Spearman rank-order correlation coefficient: 0.682, p < 0.001). The same has been observed when considering 10-km buffer areas for comparative purposes. In other words, population densities calculated at the watershed scale, easily obtained in using secondary and freely available databases, satisfactorily document demographic pressures at local scales, independently of the basins. Population densities calculated at the watershed scale will further be used as a proxy of the global human-pressures exerted on water masses.
Possible correlations between population densities at the watershed scale and LULC at both local and watershed scales were then studied. Results indicate redundant relationships between population densities at the watershed scale and several LULC classes at both geographical scales (Table 3). Firstly, population densities appeared strongly and positively correlated with the extent of artificial areas at the watershed scale (p < 0.001), indicating that population increase induces infrastructures creation, but destroys natural environment (i.e., dense vegetation areas) as indicated by the negative coefficient (p = 0.002). All the LULC classes identified at the watershed scale, except for water surfaces, were strongly correlated with their corresponding class at local scales (p < 0.001). We have already reported the temporal persistence of global trends, which appear thus widely shared spatially, regardless of the size and location of the watersheds and SRs concerned. Additionally, a strong negative correlation existed between water surfaces at the watershed scale and the extent of dense vegetation areas at the local scale (5-km buffer area: p < 0.001). This 'inter-scales' correlation indirectly corroborates our earlier assertion that water surfaces at a regional scale are intimately associated with agriculture-prone areas. High negative correlations between 'artificial areas' and the 'dense vegetation' classes at the local scale, on the one hand, and population densities at watershed scale on another hand (first column), suggest that these two LULC classes cannot provide complementary information (in addition to population densities at watershed scale) to characterize pressures on SRs. Table 3 further reveals that cropland areas at the local scale were not correlated with population densities at the watershed scale (p = 0.927). We have seen ( Table 2) that cropland areas at the local scale were systematically and negatively correlated with sparse vegetation and, at the same time, never correlated with dense vegetation surfaces. Consequently, the cropland class determined at buffer area scales constituted the most discriminating class when considering the pool of the 37 sites of interest ( Figure 4). Owing to the absence of correlation with population densities calculated at the watershed scale, the cropland class identified at the 5-km buffer areas scale was selected as a proxy of the local human-pressures exerted on water masses. indirectly corroborates our earlier assertion that water surfaces at a regional scale are intimately associated with agriculture-prone areas. High negative correlations between 'artificial areas' and the 'dense vegetation' classes at the local scale, on the one hand, and population densities at watershed scale on another hand (first column), suggest that these two LULC classes cannot provide complementary information (in addition to population densities at watershed scale) to characterize pressures on SRs. Table 3 further reveals that cropland areas at the local scale were not correlated with population densities at the watershed scale (p = 0.927). We have seen ( Table 2) that cropland areas at the local scale were systematically and negatively correlated with sparse vegetation and, at the same time, never correlated with dense vegetation surfaces. Consequently, the cropland class determined at buffer area scales constituted the most discriminating class when considering the pool of the 37 sites of interest ( Figure 4). Owing to the absence of correlation with population densities calculated at the watershed scale, the cropland class identified at the 5-km buffer areas scale was selected as a proxy of the local human-pressures exerted on water masses.  Ultimately, the above observations led to the selection of two proxies for introducing anthropogenic pressures in the analyses and interpretations of water quality data: Ultimately, the above observations led to the selection of two proxies for introducing anthropogenic pressures in the analyses and interpretations of water quality data:

1.
Contribution of the cropland class to LULC calculated in buffer areas of 5-km radius around reservoirs, associated with local anthropogenic pressures. This is possible because synoptic information is available at this scale: (i) such a buffer area corresponds to a 78.5 km 2 area; (ii) remote sensing information is based on Landsat ETM+/OLI-based description of land use at the national scale in using 30 m × 30 m pixels; (iii) there is thus enough information for defining classes and then calculating cropland contribution.

2.
Population densities calculated at the watershed scale, associated with global anthropogenic pressures. This is possible because (i) this corresponds to the least disaggregated data level that is not an administrative division; and (ii) these data are not correlated to the data quantifying the cropland class contribution to LULC at local scales. Table 1 summarizes the main characteristics of SRs studied in 2014. Depth was <3 m for more than 85% of the studied sites, with 1.45 m as median value. Transparency was extremely contrasting, varying between 1.5 cm in Koubri and 137 cm in Bala (median value: 20 cm) Water temperature was uniformly elevated (mean temperature: 28.9 • C; coefficient of variation < 5%), and conductivity varied between 73.8 µS·cm −1 in Dissin and 480.9 µS·cm −1 in Pouytenga (median value: 115.4 µS·cm −1 ). None of the sampled SRs was stratified (i.e., temperature and conductivity values were the same across the water column), and sub-surface samples were thus representative of the whole water column. The pH was near neutral for all water samples (median value: 7.6) Spectrofluorometric assessments of phytoplankton biomass were made for 28 sites: biomass concentrations varied between 1.2 and 66.3 µg·Chl a·L −1 in Petit Bagré and Saaba respectively (median value: 20.9 µg·Chl a·L −1 ). Correlations analysis between SRs characteristics (Table 4) revealed positive significant links between the depth (measured during the campaign), the volumetric capacity of the water bodies and the surface area of their watersheds, which were correlated one to the other as well. Neither the population density calculated at watershed level (Inh·km −2 ), nor the contribution of cropland to the LULC within 5-km radius buffer areas around reservoirs (%) were correlated to the other descriptive variables (depth, environmental parameters, capacity and watershed area: p > 0.01). Observed SPM loads were extremely contrasting and varied between 5 mg·L −1 (Bala) and >3 g·L −1 (Pouytenga). SPM concentrations were fully disconnected from other environmental parameters, which were independent of each other as well. Secchi disk disappearance depth was directly controlled by SPM loads. The PIM fraction contributed predominantly to the suspended particulate loading (>50% for 82% of the sites) and appeared intimately linked to the SPM load (correlation coefficient 0.846; Table 4). Very significant inverse correlations linked SPM ( Figure 5) and PIM (not shown) to the depth of the sampled sites. These correlations are strengthened when the extremes, i.e., Bagré Lake (deepest, 11.9 m, but extremely turbid site) and Bala Lakes (the clearest despite its relative shallowness, 1.6 m), were excluded from the analysis. The observed inverse correlations between SPM values, the size of the watersheds and the capacity of sampled water masses corresponded in reality to co-variations between the two later (descriptive variables), which were both closely correlated with the depth of the sampling stations, as previously mentioned.

Water Quality Characteristics
Water 2020, 12,1967 11 of 25 conductivity at 25 °C; Chl a: sub-surface chlorophyll a concentration; Cropland: cropland class contribution within 5-km radius buffer areas around reservoirs; Surface: surface of watersheds; Pop. Density: population density at the watershed scale; Capacity: volume of studied water masses). Statistically significant results are in bold.
Observed SPM loads were extremely contrasting and varied between 5 mg·L −1 (Bala) and >3 g·L −1 (Pouytenga). SPM concentrations were fully disconnected from other environmental parameters, which were independent of each other as well. Secchi disk disappearance depth was directly controlled by SPM loads. The PIM fraction contributed predominantly to the suspended particulate loading (>50% for 82% of the sites) and appeared intimately linked to the SPM load (correlation coefficient 0.846; Table 4). Very significant inverse correlations linked SPM ( Figure 5) and PIM (not shown) to the depth of the sampled sites. These correlations are strengthened when the extremes, i.e., Bagré Lake (deepest, 11.9 m, but extremely turbid site) and Bala Lakes (the clearest despite its relative shallowness, 1.6 m), were excluded from the analysis. The observed inverse correlations between SPM values, the size of the watersheds and the capacity of sampled water masses corresponded in reality to co-variations between the two later (descriptive variables), which were both closely correlated with the depth of the sampling stations, as previously mentioned.

Impact of Anthropogenic Factors on Water Quality
Possible correlations between SRs characteristics (SPM, PIM, depth, environmental parameters, surface area, capacity) and the two selected indicators of anthropogenic pressures have been explored (Table 4). There existed a statistically very significant positive correlation between SPM values and watershed population density (correlation coefficient 0.534; Table 4), which indicated that anthropogenic pressures exerted at this scale do stimulate the turbidity of downstream water masses (i.e., high SPM loads). Poor waste management associated with high population densities, increased infrastructural development, erosive agricultural practices and subsequent expansion of bare areas are example activities that could be the cause of this. Population densities were indeed intimately but inversely correlated to the LULC class associated to dense vegetation, and positively correlated to the LULC class associated to artificial areas, in both cases at the watershed level (correlation coefficients −0.465 and 0.601, respectively; Table 3). The contribution of crop areas within 5-km radius buffer zones around reservoirs was not linked to SPM and PIM values (p = 0.157 and 0.191, respectively; Table 4). This would suggest that agricultural practices in close proximity to SRs do not directly control SPM and PIM loadings of water bodies. Neither the population density at watershed level

Impact of Anthropogenic Factors on Water Quality
Possible correlations between SRs characteristics (SPM, PIM, depth, environmental parameters, surface area, capacity) and the two selected indicators of anthropogenic pressures have been explored (Table 4). There existed a statistically very significant positive correlation between SPM values and watershed population density (correlation coefficient 0.534; Table 4), which indicated that anthropogenic pressures exerted at this scale do stimulate the turbidity of downstream water masses (i.e., high SPM loads). Poor waste management associated with high population densities, increased infrastructural development, erosive agricultural practices and subsequent expansion of bare areas are example activities that could be the cause of this. Population densities were indeed intimately but inversely correlated to the LULC class associated to dense vegetation, and positively correlated to the LULC class associated to artificial areas, in both cases at the watershed level (correlation coefficients −0.465 and 0.601, respectively; Table 3). The contribution of crop areas within 5-km radius buffer zones around reservoirs was not linked to SPM and PIM values (p = 0.157 and 0.191, respectively; Table 4). This would suggest that agricultural practices in close proximity to SRs do not directly control SPM and PIM loadings of water bodies. Neither the population density at watershed level (Inh·km −2 ) nor cropland areas within 5-km radius around reservoirs (%) were correlated to the other descriptive variables (depth, environmental parameters, capacity, surface: p > 0.01). There existed a very significant positive correlation between PIM values and, again, the population density at the watershed level (correlation coefficient 0.523; Table 4). Conversely, PIM values were not correlated with watershed's surfaces (p = 0.068; Table 4), while these surfaces were strongly associated with both SPM values and depth values. The lack of correlation between PIM values and watershed's surfaces on the one hand, and the existence of a very significant negative correlation between PIM values and depth (which is nevertheless strongly positively correlated with the watershed surface), on the other hand, underlined the role that vertical water motions and sediment resuspension might play within such shallow and well-mixed aquatic ecosystems.

Diachronic Perspective
SPM and PIM values obtained during the field trip of 2014 were compared to equivalent data (same season, same field and lab procedures) gathered about ten years earlier in 21 SRs (Figure 1). The samples were collected exactly at the same stations (GPS validation), during the same hydrological period (the heart of the dry season). Depths of the SRs measured at both dates were not statistically different (Mann-Whitney non-parametric rank-sum test, p = 0.347), and the two series have been considered as comparable. The pairs of SPM and PIM series obtained from both campaigns were not statistically different (Mann-Whitney non-parametric rank-sum test, p > 0.1) despite an apparent global increase in both parameters ( Figure 6). SPM values increased for 17 lakes out of 21, while PIM values increased for 15 lakes out of 21 ( Figure S1 Supplementary Materials).
Water 2020, 12,1967 12 of 25 (Inh·km −2 ) nor cropland areas within 5-km radius around reservoirs (%) were correlated to the other descriptive variables (depth, environmental parameters, capacity, surface: p > 0.01). There existed a very significant positive correlation between PIM values and, again, the population density at the watershed level (correlation coefficient 0.523; Table 4). Conversely, PIM values were not correlated with watershed's surfaces (p = 0.068; Table 4), while these surfaces were strongly associated with both SPM values and depth values. The lack of correlation between PIM values and watershed's surfaces on the one hand, and the existence of a very significant negative correlation between PIM values and depth (which is nevertheless strongly positively correlated with the watershed surface), on the other hand, underlined the role that vertical water motions and sediment resuspension might play within such shallow and well-mixed aquatic ecosystems.

Diachronic Perspective
SPM and PIM values obtained during the field trip of 2014 were compared to equivalent data (same season, same field and lab procedures) gathered about ten years earlier in 21 SRs (Figure 1). The samples were collected exactly at the same stations (GPS validation), during the same hydrological period (the heart of the dry season). Depths of the SRs measured at both dates were not statistically different (Mann-Whitney non-parametric rank-sum test, p = 0.347), and the two series have been considered as comparable. The pairs of SPM and PIM series obtained from both campaigns were not statistically different (Mann-Whitney non-parametric rank-sum test, p > 0.1) despite an apparent global increase in both parameters ( Figure 6  We then studied possible correlations between the rates of change between the two dates (expressed in percent of 2002's values) of SPM and PIM values and of LULC classes ( Table 5). The evolution of SPM was correlated, though insignificantly (p > 0.01), to the evolution of the water and artificial areas classes which were strongly negatively correlated one to the other (p = 0.003), while that of PIM was significantly positively correlated with the evolution of cropland areas at local scales (p = 0.006). At the same time, the evolution of dense vegetation at the local scale had a statistically significant negative correlation with the evolution of the class of artificial areas (p = 0.002). This We then studied possible correlations between the rates of change between the two dates (expressed in percent of 2002's values) of SPM and PIM values and of LULC classes ( Table 5). The evolution of SPM was correlated, though insignificantly (p > 0.01), to the evolution of the water and artificial areas classes which were strongly negatively correlated one to the other (p = 0.003), while that of PIM was significantly positively correlated with the evolution of cropland areas at local scales (p = 0.006). At the same time, the evolution of dense vegetation at the local scale had a statistically significant negative correlation with the evolution of the class of artificial areas (p = 0.002). This suggests that, in the period between the two campaigns, the contribution of the dense vegetation class to LULC within 5-km buffer of SRs has been significantly eroded, and this benefited the contribution of the class of artificial areas, which increased. The evolution of the class of artificial areas is believed to have marginally stimulated the increase in SPM load between the two periods (p = 0.037). However, the augmentation of the PIM contribution to the SPM load during the same period might firstly be associated with the expansion of cropland surfaces in the immediate vicinity of SRs (Figure 8), which could eventually result in an increase in artificial areas. This can be related to poor agricultural practices that increase land degradation and promote erosion sensitivity.
Water 2020, 12, 1967 13 of 25 suggests that, in the period between the two campaigns, the contribution of the dense vegetation class to LULC within 5-km buffer of SRs has been significantly eroded, and this benefited the contribution of the class of artificial areas, which increased. The evolution of the class of artificial areas is believed to have marginally stimulated the increase in SPM load between the two periods (p = 0.037). However, the augmentation of the PIM contribution to the SPM load during the same period might firstly be associated with the expansion of cropland surfaces in the immediate vicinity of SRs (Figure 8), which could eventually result in an increase in artificial areas. This can be related to poor agricultural practices that increase land degradation and promote erosion sensitivity.

LULC Changes
Landsat analyses performed for the years 2002 and 2014 unambiguously demonstrated an increase of human induced pressures at different spatial scales around SRs in Burkina Faso, with an increase in cropland and artificial areas and a general decrease of natural spaces. The conversion of savanna to cropland in this region has already been documented by several authors [14][15][16]70,71]. Farming remains characterized by low inputs, and production increase relies most often on field expansion rather than on productivity strengthening. Such a practice, combined with population growth, explain the conversion of natural vegetation to cropland area, as observed in the present study in Burkina Faso, and by several authors in the neighboring Northern Ghana [72][73][74]. Important inter-provincial and international migrations [75] exacerbate demographic dynamics and then associated LULC changes in migration destination areas [76,77], as this is the case in the southern and central regions of Burkina Faso. In addition, fallowing practice has reduced throughout the region [78] resulting in soil quality decline and fertility loss, which favors the augmentation of bare surfaces [15]. Expansion of human settlement across the country also contributed to the extent of artificial areas and bare surface at the expense of natural vegetation, as observed elsewhere [79]. Other anthropogenic activities, such as wood harvesting for domestic uses or charcoal production, may contribute to the observed landscape degradation. Ultimately, the diminution of natural vegetation may be also exacerbated by climate change, which acts synergistically with anthropogenic pressures [80][81][82].

LULC Changes
Landsat analyses performed for the years 2002 and 2014 unambiguously demonstrated an increase of human induced pressures at different spatial scales around SRs in Burkina Faso, with an increase in cropland and artificial areas and a general decrease of natural spaces. The conversion of savanna to cropland in this region has already been documented by several authors [14][15][16]70,71]. Farming remains characterized by low inputs, and production increase relies most often on field expansion rather than on productivity strengthening. Such a practice, combined with population growth, explain the conversion of natural vegetation to cropland area, as observed in the present study in Burkina Faso, and by several authors in the neighboring Northern Ghana [72][73][74]. Important inter-provincial and international migrations [75] exacerbate demographic dynamics and then associated LULC changes in migration destination areas [76,77], as this is the case in the southern and central regions of Burkina Faso. In addition, fallowing practice has reduced throughout the region [78] resulting in soil quality decline and fertility loss, which favors the augmentation of bare surfaces [15]. Expansion of human settlement across the country also contributed to the extent of artificial areas and bare surface at the expense of natural vegetation, as observed elsewhere [79]. Other anthropogenic activities, such as wood harvesting for domestic uses or charcoal production, may contribute to the observed landscape degradation. Ultimately, the diminution of natural vegetation may be also exacerbated by climate change, which acts synergistically with anthropogenic pressures [80][81][82].  [83]. They were substantially more elevated than values observed in Burkina Faso during the dry season of 2004 in a set of 10 (large) reservoirs (median value: 16.0 mg·L −1 ; [84]). These values were also higher than SPM values observed in the lentic water bodies of the Inner Niger Delta in Mali during a two years survey, which varied between 7.4 and 34.5 mg·L −1 depending on the season and on the location [85]. Beyond median values, large variations in SPM concentrations have been encountered (three orders of magnitude: from some mg·L −1 in Bala to more than 3 g·L −1 in Pouytenga). Undoubtedly, local conditions including the location (a natural reserve and an urban site, respectively) and the mean depth of the SR (160 cm and 60 cm, respectively) did influence the observed values. Meanwhile, SPM concentrations appeared fully disconnected from other environmental parameters.

SPM Concentration Measurements
Turbidity assessment provides relevant information related to the ecological functioning of water masses and their interactions with the environment. SPM and PIM concentrations within a water body depend on both long-term (interannual and seasonal) trends (e.g., runoff) but also on short-term events (meteorological, human or animal-induced disturbances), and may thus vary in time and space [86]. Here, SPM assessments were performed on samples collected during the heart of the dry season, i.e., long weeks after the last runoffs and before the onset of the next rainy season. During the dry season, particulate matter transported by floods do settle and water transparency increases slowly in SRs, while the water level decreases progressively [87]. At this time of the year, SPM concentrations reach their annual minimum (see Figure 5 in [17]). SPM concentrations measured during our field campaigns corresponded thus to a sort of baseline, which integrated both the external inputs and the internal processes of production [88]. None of the studied reservoirs were stratified, as previously noted by other studies in this area [35]. SPM and PIM values, associated to sub-surface grab samples, exemplified thus the intrinsic properties of the whole water body.
SPM samples were systematically collected early in the morning (median value 07:40 am in 2014). This timing does not correspond to a period of intense frequentation of SRs by humans for domestic uses or by livestock herds for drinking: we assume that SRs had not been artificially disturbed before our sample collection. On another note, it is known that during the dry season, gusts of very violent stormy winds (i.e., reaching >30 m·s −1 , [89]) occur stochastically, inducing rapid and total mixing of the water column and immediate pulses of resuspended sediments. Such events were not encountered during our sampling collection. Ultimately, there exists a well-known diurnal pattern of northeast trade winds during the dry season, with wind blowing vigorously during the hot hours of the day, but never at night and early morning. Intense wind-induced sediment resuspension linked to such pattern were not likely to occur owing to our sampling timeframe. However, the duration of the impact of such regular and repetitive energy input into the water column may largely exceed hourly patterns, and one may not discard the hypothesis of a turbid and inorganic-dominated environment in SRs largely controlled by daily wind-induced motions [90,91]. Moreover, shallow sites that are intensively frequented never cleared up and always exhibit elevated turbidity, firstly because of the earthy or sandy nature of their soils. This explains why PIM contribution to SPM load was dominant compared to organic matter, and increased when SPM concentration increased. Such a close relationship is a classical observation in shallow water bodies [92,93], where it is most often associated with the remobilization of already settled inorganic particles whatever their origins [91,94]. Repetitive sediment resuspensions [95] induce random sedimentation over large areas [96]. The main consequence is that in absence of progressive accumulation to deeper sites, the same fine material remains liable to resuspension and contributes to the maintenance of elevated turbidity. Finally, [95] underlined the importance of sediment transfer from the edge of lakes, firstly attributed to the wave-induced shore erosion and redistribution of these coastal sediments. Such process might be of particular importance in SRs, owing to the increasing use of their immediate margins for vegetable production that hugely disturb the physical stability of shoreline soils. Gardening has indeed often been made responsible for SRs siltation in Burkina Faso, even if this has never been scientifically demonstrated. Meanwhile this activity is severely denounced as unsustainable agricultural practice in terms of soil management and erosion control [97].

Synchronic Analysis
We observed that SPM loads measured in 2014 were intimately associated to population densities calculated at the watershed scale for this same year. Contamination of surface waters by human activities follows two ways: (1) by point sources (e.g., sewage treatment discharge and storm-water runoff), and (2) by non-point sources or diffuse contaminations (e.g., runoff from agricultural areas) that are often difficult to detect because they generally encompass large areas in drainage basins [98]. Wide contributing areas are thus needed for being included in analyses, in order to take into account distant sources of possible disturbances [99]. Independent of latitude and socio-economical contexts, the watershed scale is often indicated as a key spatial scale on which land use affects water quality parameters, especially turbidity and nutrients [100][101][102]. We have adopted population density at watershed level as a proxy of overall pressures exerted by human populations on their environment, owing to the huge autocorrelations observed between population densities and LULC' classes distribution at this same geographical scale. Indeed, increasing population, urban expansion and agricultural intensification exacerbate soil degradation, including erosion and sediment transports, in sub-Saharan Africa [103]. Clearing of vegetation leaves land bare between cultivation cycles that accelerate soil erosion; the replacement of forest by cropland is historically used as a landscape scale indicator of degradation [68]. It is well-known that the conversion of natural landscapes into cropland promotes their sensitivity to water and wind erosion [79,104], thus indirectly contributing to SRs turbidity and siltation. Previous studies have indeed recorded increased levels of SPM in agricultural land use areas [105]. Elevated rates of soil loss, associated with poor management of cropping (intensification) and grazing (overexploitation of natural resources) contribute to downstream sedimentation and degradation of water bodies [53,106]. Moreover, land use and land use change are known to alter particle size characteristics of SPM in aquatic systems [107]. Intense land use practices typically produce greater quantities of fine (<5 µm) inorganic particles relative to other land uses [108]. Their positive-buoyancy obviously contributes to increased turbidity.

Diachronic Perspective
Twenty-one of the 37 SRs studied in 2014 had been already sampled ten years before, allowing a diachronic perspective on their properties. We found that the changes of the PIM contribution to the SPM loadings between the two dates were sensitively associated to the evolution of cropland within 5-km buffer areas within the same period. Indeed, the intrinsic properties of watersheds and buffer areas around SRs (geology, slope, etc.) did not change between the two dates, in contrast to land use imposed by demographic dynamics. Similarly, [17] already mentioned an increase in turbidity of the Bagré Lake between 2000 and 2015, firstly attributed by the authors to 'an increase of bare soils and croplands areas in this region' as observed by our study. Landscape changes (LULC, demography) are known to induce changes in the ecology and health status of aquatic ecosystems [36,37,109]. Studies found that higher percentages of urban or agricultural lands are usually associated with poor water quality, while higher percentages of natural vegetation and notably forest are related to better water quality [110]. Controls exerted by landscape changes on water quality have been classically documented, and numerous studies have highlighted in particular the role of cropland on sediment loads within water masses [111,112]. For example, [113] demonstrated the existence of a close link between the proportion of cropland and the percentage of inorganic matter within the sediments of 22 natural marshes in Ontario. Concentration of clays systematically increases with denaturation of landscape and its conversion towards cropland [105,114], then justifying the diachronic association between enhanced PIM concentrations and increased cropland areas as observed here.

Defining Scales of Interest
The issue of whether land use near water bodies (local scale) is a better predictor of water quality than land use over the entire watershed has remained unsettled for a long time [115][116][117]. Two spatial scales (including the whole watershed and 1 to 5 km riparian buffers) are often used to interpret the relationship between the characteristics of watershed and water quality measures, but which one is better is still in big arguments [118,119]. Because land uses are not randomly distributed throughout basins, changing the scale of the contributing area modifies the relative proportion of various land uses, the perceived dominant land use, and the perceived rate of land use change [120,121]. For example, [115] reported that nutrient and fecal coliforms concentrations were better predicted by local riparian land use rather than watershed land use in small low-order streams in pasture watersheds in New Zealand. [116] also concluded that the impact of land use on water quality is stronger in the immediate buffers around water masses than that in the watersheds. These findings may be related to the differences in watershed' land use patterns or to the water quality parameters of interest. Moreover, the relationships between land use changes and water quality indicators are known to be scale dependent [101,122]. There is nevertheless a growing consensus acknowledging that management options need to address both processes at large scales, (e.g., watershed-wide forest land cover and urban expansion), and processes that develop on small scales (e.g., agricultural activities and riparian vegetation near water bodies) [101,123]. This has been for example recently highlighted for lakes and reservoirs of Northeast Brazil, a climatically comparable region of Burkina Faso [124]. Authors emphasized that: "( . . . ) to safeguard water quality, not only the protection and conservation of nearby areas of the aquatic systems are important, but also it is fundamental to consider the land use characteristics of the whole catchment area ( . . . )" and that "This is especially relevant for small and shallow lakes." ( [124] and references herein). Simpultaneously considering watersheds ('global pressures', here exemplified by population densities) and small buffer areas around water masses ('local pressures', here the contribution of cropland within 5-km buffer areas) as we did, allows taking into account both the influence of upstream inputs and of lateral 'immediate' inputs.
The number of SRs involved in the comparison of water quality between 2004/5 and 2015 (21) is admittedly a limitation of this study. Nonetheless, results obtained are promising and provides basis for stimulating further works, in including more SRs and more water quality indicators. This recommendation aims hopefully at pave the way for future research.

Conclusions
The possible impacts of anthropogenic pressures (LULC and demography) on ecological properties of Sub-Saharan small reservoirs have never been investigated. This study sought to fill this gap by relating demography and LULC changes at different geographical scales and water quality changes in a decade. We firstly found important redundancy between downloaded population data and remote-sensing deduced LULC information. This allows a cross-validation of the two sources of information, and validates the relevance of the approach. Two complementary indicators have been identified that reflect the nature and intensity of pressures exerted by humans on their environment, and, indirectly, on aquatic ecosystems. The involvement of population densities calculated at the watershed scale is the most inclusive proxy that can be defined for encompassing global human pressures in their diversity. This is also the most easily accessible secondary data (open access gridded population density data are available on the internet), which needs 'only' the delineation of watershed to be included in future analyses. This proxy can thus be used as an indicator of the trends associated to 'the overall pressure exerted by humans on the landscapes', here exemplified by the positive relationship linking SPM loads within SRs and population densities of their basins. The proportion of cropland within small buffer areas around reservoirs (5-km) is not correlated to population density. The involvement of this unique LULC can thus be considered as a complementary, easily available and valid proxy of the intensity of 'local anthropogenic pressures immediately around water masses', here illustrated by the positive relationship linking small-scale LULC evolutions in a decade, and changes in PIM contribution to SPM loads of SRs over the same period.

Recommendations
The implementation of soil and water conservation measures at the watershed scale in order to increase agricultural rainwater productivity and to limit efficiently soil erosion and sediment transport, has to be supported: beneficial for farmers, such measures do also protect downstream aquatic ecosystems [125][126][127]. Moreover, at the catchment scale, a series of evidences suggests that LULC-diversified landscapes reduce fluxes of nutrients and contaminants towards water bodies, when compared to crop-specialized landscapes, and mitigate efficiently the impacts of diffuse pollution [128][129][130]. Proactive investments for maintaining watersheds and landscapes diversity have to be promoted. Meanwhile, protection measures immediately around water masses have also to be encouraged. This point does correspond in reality to the official and legal posture in Burkina Faso (article 7 of the Law n • 014/96/ADP, in date of 23 May 1996) which indicates that there should exist a 100 m-easement around water bodies where human activities are not allowed. Protection perimeters where agriculture is excluded, is classically recommended for alleviating erosion and then siltation [97]. Forested riparian buffer strips are known for their recovery functions, in catching nutrients and sediments and then protecting water masses [131,132]. Their crucial role for the pollination of adjacent agriculture fields has also been recently demonstrated [133]. Tree planting and re-greening the immediate vicinity of SRs [134] could thus be considered as a valuable perspective for increasing the multiple ecosystem services associated with SRs, as recently experimented in Burkina Faso [135]. Increasing the vegetation within buffer zones around water bodies in order to enhance their quality has indeed been urgently pleaded [136]. Such a nature-based solution could be a proposal for the management of the immediate neighborhood of SRs. Farmers are the primary actors on SRs landscapes and their perspective and expectations have to be incorporated in management strategies, particularly when strategies for reversing degradation and improving food security are implemented [103]. It remains that there exists a fundamental nexus between water resource protection and the intensification of agricultural practices, in terms of spatial extension and in terms of inputs (e.g., fertilizers and pesticides). This obligatorily requires in a near future the implementation of explicit conservation policies focused on SRs, in order to maintain the many ecosystem services they do provide to rural populations [137]. Research on the impact of management practices on erosion control, namely regarding agricultural practices, reforestation and protection of riparian areas, is also urgently needed [138].