Land Cover Change and Water Quality: How Remote Sensing Can Help Understand Driver–Impact Relations in the Lake Titicaca Basin

Abstract

The increase of human interventions and developments are modifying the land use/land cover (LULC) of the global landscape, thus severely affecting the water quality of rivers and lakes. Appropriate management and effective policy developments are required to deal with the problems of surface water contamination around the globe. However, spatiotemporal variations of water quality and its complex relation with land cover (LC) changes, challenge adequate water resources management. In this study, we explored the use of remote sensing to relate LC change in the Katari River Basin (KRB) located in the Bolivian Andes and water quality on the shores of Lake Titicaca, in order to support water management. An unsupervised classification of Landsat 7 satellite images and trajectory analysis was applied to understand the modifications of LC through time. In addition, water-quality indicators at the outlet of the basin were retrieved from remote-sensing images and its temporal behavior was analyzed. The results show that the expansion of urban areas is the predominant environmental driver in the KRB, which has great impact on the water quality of Lake Titicaca. We conclude that there is a strong link between the rapid growth of urban and industrial areas with the detriment of river and lake water quality. This case study shows how remote sensing can help understand driver–impact relations.

1. Introduction

Urban population growth, increase in agriculture, mining activities, industrial developments, among others, influence the environment, and are of global concern because of their socio-ecological impacts [1]. Large land use modifications and management practices are key influencing factors altering the normal functioning of the hydrological system, causing alterations in runoff and water quality [2,3,4,5].

Land cover (LC) mapping is one of the major important applications of satellite-based remote-sensing data as it allows to study the LC dynamics with a wide range of spatial and temporal coverage for a variety of resource-based applications [6,7,8]. The method employed for the LC mapping for the research reported on in this article was unsupervised classification. This method is based on a chosen algorithm that finds a specified number of statistical clusters in a multispectral space, and it is usually employed when there is not enough data for training or when there is no prior knowledge of the cover in the study site [6,9]. LC change is important for monitoring the local resources and environment, as the historical data makes it possible to analyze spatiotemporal patterns of environmental elements and the impacts of human activities [10]. Tracking the changes of LC throughout the years can help to understand the development of the region and provide references for resource management [9]. Land use and land cover (LULC) change studies are imperative to provide sustainable land resources use, rehabilitation measures, and evidence-based support to improve management [11]. Furthermore, LULC changes greatly influence hydrologic processes and fluxes (i.e., surface runoff and stream flows) due to variations in the physical characteristics of the land surface, soil, and vegetation [11,12,13]. In some cases, hydrological models (i.e., VIC, SWAT, HEC-HMS, DHSVM, InVEST, among others) have been employed to determine the effects of LULC changes in hydrology, in which an increase in surface runoff and a decrease in evapotranspiration and infiltration are the most common responses [12,14,15,16,17]. Trajectory analysis is a method used for LC change research based on the time series of each pixel, which is used to detect trends [18] and quantify land cover transformations through time. Consistent mapping of LULC and its changes over time is necessary to measure the impact of local actions and global indicators such as the SDGs [19].

Satellite-based remote sensing also has the potential to provide spatiotemporal information regarding water-quality conditions over large areas [20,21,22]. In this research, remotely sensed data have been used in combination with algorithms to assess water quality, allowing, e.g., chlorophyll-a and turbidity mapping in surface waters [22,23,24,25,26]. These retrievals have been successfully used for monitoring inland and coastal waters [27], and to even measure changes in water quality during the recent lockdown due to the pandemic [28]. Remotely sensed chlorophyll-a retrieval can provide complementary information to in situ data to support Water Framework Directive monitoring requirements [29]. Turbidity (TUR) is a measure of the degree to which the water loses its transparency due to the presence of suspended particulates [30]. Whereas, chlorophyll-a (CHL) is a green pigment found in plants and is considered an indicator of phytoplankton abundance and biomass in surface waters [20]. Both parameters are indicators of the overall quality of surface waters. Furthermore, remote-sensing data is also used for the identification, analysis and mapping of vegetation dynamics based on vegetation indices [7,31]. These indices usually consist of a combination of two-, three-, or four-band formulae [31,32,33], and in this research, they will be used as a proxy for measuring the eutrophication in Lake Titicaca.

Changes in LULC have the potential to affect surface water quality at diverse scales [4]. LULC change analysis at basin and sub-basin scales is rare, which is unfortunate considering its usefulness and importance for the design of appropriate land management practices, strategies, and policies that best fit local conditions [34]. A study on the link between LULC and water quality in river basins has the potential to identify pollution that can provide stakeholders with information regarding potentially harmful effects [21].

The implementation of effective water-quality management strategies requires an in-depth understanding of the activities that are causing such environmental modifications [35,36]. This understanding can only be achieved by having regular, affordable, and non-invasive approaches for assessing water quality, with adequate spatial and temporal coverage [22]. Thus, to adopt effective actions for the protection of water bodies, an analysis of the cause–effect relation between LULC and water quality requires special attention [37].

River basins constitute complex socio-ecological systems with different levels of interrelations between their social, economic, physical, chemical, and biological components [38]. On one hand, LC changes are very site-specific, based on natural and anthropogenic developments, but on the other hand, the effects of such modifications tend to manifest in a different local spatial scale [3,21,37,38]. Therefore, in land and water management it is important to take off-site effects of changes within the basin (on-site) processes into account when deciding on a specific intervention [35]. Here, the concept of “connectivity” illustrates that although management decisions are made within a previously defined area (namely the basin), these decisions have effects both within and outside this geographic unit [35]. Therefore, this connectivity must be valued when making management decisions, to ensure the sustainability of the landscape [35,36,38]. Effective policymaking requires a detailed understanding of the dynamics and cause–effect relations at each site [4,34,35]. Here lies the novelty of our research, which is based on the principle of understanding the relationship between changes in land cover (on-site) and their effect on the quality of the waters of Lake Titicaca (off-site). This is to detect and highlight areas within the basin where measures can be employed that have an effect not only within the river basin, but also in the lake. This will provide the public authorities and decision-makers with the tools and knowledge to implement more sustainable LULC planning and management.

This study assesses the spatial and temporal dynamics of LC change in the KRB in Bolivia, South America. This basin incorporates the presence of mining, urban, industrial, and agricultural developments [39,40]. These human developments largely modified the surface water quality of the system with severe consequences for local indigenous communities situated in the downstream region [38,41]. Furthermore, this research explores driver–impact relations with the water quality of Lake Titicaca, which represents the most important water resource in the Andes region [38,41].

The traditional approach to solving water (quality) problems are technological fixes and end-of-pipe solutions [42], which is the case for most of the policy responses in the KRB that only tackle the “symptoms” of the problems [38]. To make the transition to integrated water management, there has to be a shift towards a more thoughtful attitude that involves developing integrated problem-solving [42]. Thus, this research represents a baseline to identify areas prone to water-quality problems and allocate resources accordingly as increasing levels of pollution call for the active management of the lake and its basin [35,43].

This study, using the KRB as the study area, evaluated and analyzed the spatiotemporal variations of land cover from 2006 to 2018 and the impacts of LULC changes on the water quality of Lake Titicaca. Furthermore, this study hypothesized that remote sensing can be used to attribute driver–impact relations of lake water quality. The research outline was as followed: (a) quantifying the land cover changes in the basin; (b) detecting and quantifying eutrophication extent in Lake Titicaca with vegetation indices as proxy; and (c) calculating trends in the spatiotemporal changes that were identified. The innovation of this study lies in the use of satellite-based remote sensing to relate LC changes and lake water quality. The results from this research may provide scientific reference for quantitative assessment, effective management, and sustainable development of water resources in the KRB located in the Bolivian Andes.

2. Materials and Methods

2.1. Case Study

The KRB is located in the Andes region near La Paz, which is the capital of Bolivia [44]. The basin spreads over an area of 3006 km² with elevations ranging from 3800 to 5200 m above sea level (m.a.s.l.) [45]. The daily temperature ranges between −5 and 15 °C, which is mainly controlled by solar radiation [41]. The climatic conditions are typical of the Andean highlands. Average annual precipitation is observed to be between 470 to 742 mm, with the highest precipitation in the northeastern, mountainous, and salty areas of the basin, due to the availability of evaporated water from the lake [45]. The dry season starts from April to September and the rainy season from October to March [39,41,46].

Approximately 11% of the Bolivian population lives in the KRB, which represents only 0.27% of the Bolivian territory, making it the most populated river basin in the country [40,41]. The most important urban sectors are the cities of El Alto and Viacha [47] since these cities represent 96.4% of the total population settled within the basin. El Alto has become one of the fastest-growing cities in Latin America [41,48], as indigenous peoples have left rural areas to seek employment and livelihoods [49]. Despite this massive population growth, the necessary sanitation infrastructure developments have not been carried out. Around 45% of the urban wastewater is released to the surface water bodies crossing the city of El Alto. It is estimated that El Alto’s discharges represent seven times the natural river flow [41].

The basin also has more than 5000 industrial facilities situated in the urban areas of El Alto and Viacha, most of them with no wastewater treatment plants [41]. Additionally, in the highest part of the basin, there are contamination problems due to the large presence of mining waste [40,41,50]. On the other hand, in the lower parts of the basin, agricultural and livestock activities are carried out and their residues are carriers of organic and microbial pollution that negatively affects human and animal health [40].

The KRB hydrology consists of 4 main rivers that flow through the basin: the Seco, Seke, Pallina, and Katari Rivers [41]. The Seco and Seke rivers originate at the Cordillera Real. The upstream section of the Seke river discharges into Milluni Lake, which is affected by mining drainage [38,39,40,41]. Downstream, both rivers flow through El Alto, where they receive ineffectively treated and untreated wastewaters from households and industries. Later, these rivers discharge into the Pallina river, which flows through the city of Viacha. Downstream of the city, the Pallina river flows into the Katari river, which crosses its way through agricultural areas. The Katari river then goes through Cohana Bay and finally discharges into Lake Titicaca, with an average river basin outflow at the discharge point of 7.7 m³/s [50]. The location and course of the aforementioned rivers can be seen in Figure 1, as well as the several municipalities that share jurisdiction over the basin. The discharge of some of the rivers in the basin experience rises by a factor of ~7 from dry to wet season, representing the high seasonal fluctuations [40].

Lake Titicaca is the second largest lake in South America, the highest lake in the world, and the most important water resource in the Andes Region [38,41]. Apart from being an essential water resource for the region, Lake Titicaca has a retention time of approximately 1000 years [51]. High levels of contamination and eutrophication have been detected in the lake, affecting the local population and causing algal blooms [41,43,47]. Therefore, the quality of the water that flows into it must be carefully monitored and managed.

The basin first caught the attention of the Bolivian national authorities when the arrival of a large mass of contaminated water to the lake was reported in 2001 [40]. This event modified the color of the water and its turbidity as a result of heavy rains and resulted in the death of 2 tons of birds, fish, and frogs [40]. As a result, for more than 15 years, local, regional, national, and international agencies have invested resources in the creation and implementation of policies to improve the situation of the basin [38]. This basin is a suitable case study for the development of this research in which paradigms shifts are needed in the water management policies as previous end-of-pipe approaches have not achieved the desired results.

2.2. Land Cover Change Analysis

Landsat 7 images were used to generate land cover maps from 2006 to 2018, using a three-year time step. Subsequently, these maps were used to perform a trajectory analysis to detect the LC change dynamics in the basin. The Mann–Kendall trend test and Sen’s slopes estimator were then applied to the outcomes.

2.2.1. Basin Area and Satellite Images

The basin area was obtained using Global Mapper v18.0 and the 30-m Digital Elevation Model (DEM) downloaded from GeoBolivia. Sub-basins were generated with a minimum stream drainage area of 400 km² and with a flow tracing at the outlet of the basin. In addition, stream lengths of less than 500 m were discarded, and a depression filling of 50 m was applied to facilitate the flow network creation.

Due to high seasonal contrast in the rivers of the basin [50], an LC change assessment every 3 years was performed for the dry and wet seasons. Level-1 satellite images were obtained from the United States Geological Survey website, the Landsat 7 satellite was chosen as it had been in operation during the entire study period. Satellite images corresponding to Path 1, Rows 71, and 72 were used due to the geographical position of the basin, with a spatial resolution of 30 m. A total of 20 satellite images were chosen for this study, and priority was given to images with a cloud cover of less than 10% (a table detailing the satellite imagery used can be found in the Supplementary Materials). This significantly reduced the availability of images and made it difficult to analyze similar dates. Despite these efforts, certain areas with cloud cover and shadows had to be masked in the wet season maps, thus, the accuracy assessment and trajectory analysis were only performed using the dry season maps. In Table 1, a summary of all the downloaded data used in this research can be found.

QGIS v. 3.16.7 was used for the preprocessing, which consisted of (a) correction of the Scan-Line error, (b) clipping the images to the desired extent, and (c) atmospheric correction and conversion to reflectance. The Scan-Line error was corrected using the “Gap Mask” files for each image and the “Fill no data” feature in QGIS. The clipping of the images, the DOS1 atmospheric correction, and the conversion from Digital Numbers (DN) to reflectance values were performed using the Semi-Automatic Classification (SCP) Plugin version 7.8.15 Matera for QGIS. According to the SCP plugin manual, the calculation of the Top of Atmosphere (TOA) reflectance was made using the equation proposed by NASA:

(1)${\rho }_p=\left(\pi \ast L_{\lambda }\ast d^2\right)/\left(ESU{N}_{\lambda }\ast cos{\theta }_s\right)$

where, L_λ = Spectral radiance at the sensor’s aperture (at-satellite radiance), d = Earth–Sun distance in astronomical units (provided in the metadata files), ESUN_λ = mean solar exo-atmospheric irradiances, and θ_s = solar zenith angle in degrees, which is equal to θ_s = 90° − θ_e, where θ_e is the Sun elevation.

Similarly, the plugin performed the atmospheric correction by using the Dark Object Substraction (DOS) following this equation:

(2)$L_p=L_{min}-L_{DO1\%}$

where, L_min = radiance obtained with that digital count value (DN_min), L_DO1% = radiance of Dark Object, and is assumed to have a reflectance value of 0.01. Atmospheric correction consists of the removal of the scattering and absorption effects from the atmosphere to obtain surface reflectance that characterizes surface properties [52]. For further details on the calculations performed in the plugin and the equations used, please refer to its manual and related literature [53].

2.2.2. LC Classification

An unsupervised classification was performed for this study due to the absence of historical data that could be used as training pixels for a supervised classification. The Iterative Self-organizing Data Analysis Technique (ISODATA) was used for clustering. This method has some refinements [54] and consists of (a) assigning arbitrary initial values to clusters, (b) classifying each pixel to the nearby cluster, (c) calculating the cluster mean of all pixels in one cluster, (d) repeating b and c steps continues until the “change” between the iteration is small, and (e) improving the ISODATA by splitting and merging clusters based on a certain threshold [55]. This technique has been applied in several cases showing satisfactory classification outcomes [31,32,56].

The clustering feature of the SCP Plugin was used for the classification and a threshold of 0.01 was defined as the distance required for merging clusters, as recommended by the plugin developer [53]. The number of classes was set as 15 to ensure a good differentiation between them and to avoid misclassification in urban areas due to similarities of spectral signatures with other classes [57].

The Normalized Difference Vegetation Index (NDVI) was calculated for all the images to assist in the visual analysis and reclassification of the land cover classes, as it is used to identify vegetation, understand vegetation density, and assess changes in plant health [58]. The Raster Calculator in QGIS was used for the calculation of the NDVI, by using the formula proposed by USGS (2021) for Landsat 7:

(3)$\mathrm{NDVI}=\left(\mathrm{Band} 4-\mathrm{Band} 3\right)/\left(\mathrm{Band} 4+\mathrm{Band} 3\right)$

Based on the NDVI, Google Earth [57], satellite images, and knowledge about the basin, the reclassification was done in the SCP Plugin. A total of six land cover classes were identified: (1) water, (2) built-up, (3) vegetation, (4) barren Land, (5) shrubs/grass, and (6) snow cover.

2.2.3. Post-Classification Enhancement

Land cover classes like water, snow cover, and built-up usually present similar spectral signatures [57], which lead to misclassification of various pixels. To correct this error, further post-classification enhancement was performed. Cloud masking was done using the SCP plugin to avoid misclassifications in the map. Based on the NDVI images, Google Earth historical imagery, and the true- and false-color images, the airport and the urban areas were delineated as suggested by Manandhar et al. [57]. These areas were then incorporated into the final classification to obtain more accurate land cover maps of the basin.

2.2.4. Accuracy Assessment

The historical imagery from Google Earth offers high-resolution satellite imagery in different dates at many places [8]. Thus, 51 polygons were created throughout the basin, and Google Earth was used to assign the class corresponding to each polygon throughout the time period.

These polygons and each classification were compared using the accuracy feature of the SCP Plugin. Only the dry season images were analyzed, as the wet season images had areas with “No Data” values due to cloud masking.

2.2.5. Trajectory Analysis

Land cover change trajectories were determined by intersecting the dry season land cover maps in QGIS. The intersection was done following a chronological order of the land cover maps, in which all the pixel histories were established, and its areas calculated.

To analyze temporal human impact on the basin, the methodology proposed by Zhou et al. [59] was used to group all the trajectories into three main categories: (a) unchanged, (b) human-induced change, and (c) natural variability.

• Unchanged category: Trajectories indicating the same land cover class found in the pixels over the 13 years. Ambiguous cases in which the pixel had a different category in just one of the years were also incorporated in this category, considering possible classification errors that could have created a once-only false class in the trajectory (unless built-up class was present as the final class in 2018, as it would indicate a recent human-induced change).

• Human-induced category: Trajectories in which the built-up class was present, indicating decisive changes due to population growth and urban expansion in the basin.

• Natural variability category: Indecisive changes due to natural processes or minor human activities such as light grazing. It usually consisted of oscillations between the snow cover and water classes due to the ice melting or floods, or between vegetation, shrubs/grass, and barren land due to the natural cycle of vegetation in the basin.

The trajectory analysis was only performed for the dry season images due to the cloud masking in the wet season images, which resulted in the loss of information that could generate inaccurate results. Once the areas for each land cover class were obtained, these data were analyzed using the Mann–Kendall test to determine trends in the data. Similarly, Sen’s slope were calculated to determine the magnitude of these trends. Both calculations were performed using the Python package PyMannKendall [60].

2.3. Water-Quality Analysis Using Remote-Sensing Data

Vegetation indices such as the NDVI are widely used across the world for monitoring vegetation cover changes due to their simple estimation and easy availability at different spatial and temporal resolutions [61]. Green biomass (algae and aquatic plants) has higher reflectance in the near-infrared waveband than in the red waveband, therefore if a water pixel has a high level of green biomass, its NDVI value will also be high [62]. Despite its advantages, NDVI is still affected by atmospheric, soil, and background vegetation cover feedback conditions [61,63].

Therefore, as an alternative, we have also utilized the Soil Adjusted Vegetation Index (SAVI), as it was established to improve the sensitivity of NDVI to soil backgrounds [64]. The equation for SAVI adapted for Landsat 7 satellite images is provided by USGS [65]:

(4)$\mathrm{SAVI}=\left(\mathrm{Band} 4-\mathrm{Band} 3)/(\mathrm{Band} 4+\mathrm{Band} 3+0.5\right)\ast 1.5$

In the case of aquatic ecosystems, the background signal is mainly due to the water component [33], therefore, as we performed these calculations for Lake Titicaca, I Normalized Difference Aquatic Vegetation Index (NDAVI) was also determined. The NDAVI is a vegetation index specifically targeted at aquatic vegetation and aims to reduce the background signal due to the water [66]. This index comes as an adaptation of the NDVI specifically targeted for aquatic vegetation [33]. The formula was adapted from Villa et al. [66] to correspond with Landsat 7 bands:

(5)$\mathrm{NDAVI}=\left(\mathrm{Band} 4-\mathrm{Band} 1)/(\mathrm{Band} 4+\mathrm{Band} 1\right)$

To provide an insight on the effects of urban growth on the water quality of Lake Titicaca, the SAVI, NDVI, and NDAVI calculations were used as a basis for identifying areas of the lake with unusually high index values for water (0.3–0.8). The SAVI, NDVI, NDAVI, and true- and false-color images were used to delineate the eutrophication patches in the lake throughout the years to calculate their extent. The Mann–Kendall test and Sen’s slope estimator were also applied to these outcomes to identify trends and magnitudes.

Similarly, the water-quality parameters of chlorophyll-a (CHL) and turbidity (TUR) available through the Terrascope platform were used as a result of the processing of satellite images obtained from Sentinel-2. The algorithms used for the processing and retrieval of these water-quality parameters are explained by De Keukelaere and Knaeps [20]. Data were only available for the last years of the study period (2015–2018) due to the launch of the Sentinel-2 space mission in mid-2015. Only one image was available for 2015, 24 for 2016, and 30 and 59 images for 2017 and 2018, respectively.

For the analysis of this information, an annual average of zonal statistics was calculated for an area of Lake Titicaca approximately 10 km from the shores (outlet of the basin). This distance was selected because it corresponds directly to Cohana Bay; more area than that analyzed for SAVI, NDVI, and NDAVI was covered to assess the state of water quality beyond this extension. In addition, it is important to mention as part of the processing of the images used, the areas considered as land coverage in the 10 m ESA WorldCover 2020 map were masked [20]. Due to the presence of aquatic vegetation and eutrophication in Cohana Bay, many areas are classified as “Herbaceous wetlands” in the World-Cover 2020 map. The algorithm used to retrieve CHL and TUR masked these regions [20], and therefore, were not available for this last analysis.

The results presented in this research are from areas of the lake even further away from the shore, thus indicating the spread of contamination. It is considered valid to assume that the concentrations of CHL in the masked areas that do not correspond to aquatic vegetation, but eutrophicated areas, are even higher.

3. Results

3.1. LC Maps and LC Change Trajectories

Figure 2 shows the resulting LC maps for the dry and wet season, in which the seasonal vegetation dynamics can be perceived, as well as the extent of the urban growth throughout time.

In Table 2, the areas for each class throughout the study period can be seen for the dry season images. By comparing the area of each class in 2018 with 2006 an increase of about 123% of built-up areas and a reduction of about 69%, 41%, and 12% could be calculated for the water, vegetation, and shrubs/grass classes, respectively. Barren land remained constant throughout the study period. However, it is important to note that all classes fluctuated between years, except for the built-up class, which was always increasing.

The overall accuracies obtained for the study period were 95.62, 91.49, 85.45, 91.58, and 94.67% respectively. All of the accuracies were above 85%, which according to Anderson et al. [67], is good enough to perform a representative trajectory analysis. The error matrix is shown in Table 3. The Mann–Kendall test showed an increasing trend (p < 0.05) for the built-up class, and the Sen’s slopes estimator determined an increasing rate of +8.023 km² per year. For the remaining land cover classes, the Mann–Kendall test showed no trend (the complete overview of results from these tests can be found in the Supplementary Materials).

A total of 1023 land cover change trajectories were identified, out of which, 65 were reclassified as human-made changes, with an area corresponding to 3.25% of the total basin. Figure 3 shows the different trajectories.

Similarly, it can be confirmed that the major changes identified in the natural variability category are those of grass/forest, with a total of 386 trajectories that account for 10% of the basin area. This trajectory is also linked to the increasing trend detected for the built-up land cover throughout the study period.

3.2. Eutrophication and Water Quality in Lake Titicaca

The results of the area that presented unusually high SAVI, NDVI, and NDAVI values (higher than 0.3) at the outlet of the KRB in Lake Titicaca can be found in Figure 4. Values are grouped by the dry and wet seasons. In most cases, it can be seen that the area with high NDAVI values exceed those calculated for NDVI. Similarly, the constant increase of the area with high NDVI and NDAVI index values over time can be noted. Additionally, the areas with high SAVI values are lower in magnitude compared to the other indices. The Mann–Kendall test showed an increasing trend (p < 0.05) for the three indices calculated for both dry and wet seasons. Similarly, Sen’s slopes estimator determined an increasing rate for all the indices (for the summarized tables, please refer to the Supplementary Materials).

First, it can be noted that the area with high NDAVI values is generally larger than the identified NDVI areas. Conversely, the areas with high SAVI values are on average 26% and 30% lower than the NDVI and NDAVI areas, respectively. Furthermore, the dry and wet season areas do not vary considerably between them. Finally, an increase of areas with high values for all indices over time is demonstrated. Comparing the areas calculated for 2018 in respect to 2006, for NDVI there was an increase of 54.1% and 45.3% for the dry and wet seasons, respectively. The same comparison for NDAVI showed an increase of 34.7% for the dry season and 32.8% for the wet season. Similarly, for the SAVI, an increase of 24.2% and 33.58% was calculated for the dry and wet seasons, respectively. In summary, an average area increase of 28.9% was estimated for the SAVI, 49.7% was obtained for NDVI, and 33.8% for NDAVI. In terms of area, for NDVI, an average growth of 2.1 km²/year was observed, 1.06 km²/year for SAVI, and 1.6 km²/year for NDAVI.

The digitized extent of the eutrophication in Lake Titicaca captured by the NDVI can be found in Figure 5. A constant behavior in the extent can be seen throughout the seasons. This indicates that the eutrophication in the outlet of the basin is so severe that not even an increase in flow due to precipitation, which causes dilution, impacts the extent (for the complete data on eutrophication extents measured by all the indices, please refer to the Supplementary Materials).

Regarding the analysis of CHL and TUR, in Figure 6, images of the distribution of concentrations in the analyzed lake extension can be observed, both correspond to measurements made on 29 September 2018. This date was selected as it presents the highest average concentrations in the studied period, 497.17 mg/m³ and 350 FNU for CHL and TUR, respectively.

The concentrations were reclassified into the ranges that are used by the Lake Titicaca Binational Autonomous Authority (LTBAA) to classify the water quality of the lake [68]. These are the OCDE trophic states for CHL and the Swedish Environmental Protection Agency levels for TUR. In Table 4, the mean concentrations of CHL and TUR for each year can be found. CHL presents a very high concentration in 2015 because only one December image was used for this year, which presented very high CHL values as it falls in the middle of the wet season. Excluding the CHL concentration in 2015, it can be seen that over the years, both average concentrations increased.

CHL experienced an increase in the concentration of 40% when comparing results from 2018 to 2016. Whereas, TUR suffered an increase in the average concentration of 155% when comparing 2018 results with that of 2016. The concentrations of 2015 were not included in these calculations as they only represent one month and not an average of the year.

4. Discussion

4.1. Trajectory Analysis of LC Changes 2006–2018

The fast urban growth mentioned in the literature [39,41,44,48] was confirmed, with an increase of 123% of built-up areas in 2018 compared to 2006 and an increasing trend. The vegetation and shrubs/grass classes in the different land cover maps showed different behavior in both dry and wet seasons from what was expected. Typically, vegetation is more abundant in the wet season due to an increase in precipitation. However, as the satellite images downloaded were taken at the beginning of each season (to ensure a cloud cover lower than 10%), the impact of precipitation, or lack thereof, could not always be captured. Nonetheless, significant changes were observed in the vegetation throughout the seasons, mainly in areas where agricultural activities are carried out. The Mann–Kendall test perceived trends only in the built-up class, which may be due to the classification deficiency found in Table 3, especially in the vegetation and shrubs/grass classes. In future research, this can be improved by employing the use of in situ data for land cover classification, improving the spatial resolution of satellite images, or by combining multiple remote-sensing data.

On the trajectory analysis, 3.25% of the watershed area was identified as having undergone changes in land cover caused by direct human intervention. Adding these changes to the urban areas present in 2006 results in an alarming urban growth of about 7.5 km² per year in the study period. The natural variability represented 10% of the basin area, which is in line with the high seasonal fluctuations in the rivers’ discharge. Some studies have investigated vegetation indices’ response to rainfall [61,69], applying this type of analysis on the KRB could provide more comprehension on seasonal fluctuations of discharges and vegetation changes.

The final trajectory map obtained can be considered of great interest for policymakers and the population in general. It is a tool that spatially displays the anthropogenic modifications in the KRB, which makes it possible to identify priority areas for the design of water management policies. This is principally for monitoring water quality in the rivers closer to the densely populated areas, because most of the streams near human settlements are the recipients of domestic and industrial wastewater discharges [70]. Therefore, the closer a river is to a populated area, the greater the possibility that it will be affected by the discharge of the aforementioned wastewater [70]. By comparing Figure 3 with Figure 5, it is evident how the change in land cover in the KRB has affected the water quality of the lake in a 12-year time period. This is the innovative aspect of this research, which not only uses remote sensing to detect changes in the environmental drivers (land cover), but also to quantify the impact they have had on aquatic resources (impacts). This type of easily visible and detectable spatial quantification has already been categorized as highly useful by the Bolivian Ministry of Environment and Water (findings were already presented to delegates from the ministry).

It is known that in this basin, a large amount of human and economic resources are being allocated to water-quality monitoring [40]; however, there are still weaknesses in public policies and their implementation [39,40,41]. This is evidenced by pollution problems identified several years ago, which are still present [41,71], and most notably in the incomplete and inconsistent water-quality database managed by the Bolivian Ministry of Environment. These problems suggest an allocation of resources without a clear delimitation of objectives [41,71] and areas of greater fragility due to the particular characteristics of LULC. The connectivity concept needs to be considered when designing and implementing future public policies on water resources to ensure the sustainability of these efforts.

4.2. Water Quality in the Outlet of the Basin

The literature suggests that water has a standard NDVI of zero or even negative values, indicating non-vegetated spaces [72]; however, near Cohana Bay, the NDVI values fluctuate between 0.3 and 0.8, which correspond to values of a densely vegetated area [73]. Likewise, low values of vegetation indices typically represent water or soil surfaces, whereas higher values represent areas of vegetation, and more specifically, in aquatic environments are indicative of algal blooms or aquatic plants [31]. When calculating the vegetation indices at the outlet of the basin, a constant extent of high index values throughout the seasons is displayed. This is further supported by the increasing trends (p < 0.05) calculated for all the vegetation indices, supporting the argument that eutrophication is expanding in the lake [41,46].

It is worth noting that Cohana Bay has a presence of aquatic vegetation [40,46], which influences these indices. Nevertheless, high levels of eutrophication have been reported in recent years from in situ measurements [40,41,46,47,71,74], therefore, the high-value areas detected by the indices were used as proxies to measure eutrophication and monitor its expansion.

As high index values are indicators of green biomass, this can be considered as an indication that the eutrophication in the outlet of the basin is so severe that not even an increase in flow due to precipitation, which causes dilution, impacts the extent. The detected growth and expansion of eutrophication coming from Cohana Bay support what has been mentioned in the literature [40].

By comparing the previously mentioned eutrophicated extensions obtained annually for NDVI and NDAVI, an average of 5% more area was measured with NDAVI than with NDVI. This finding is in line with Villa et al. [66], who revealed that NDAVI gets sensitivity results to chlorophyll content similar to NDVI’s. However, with SAVI, an average of 25–30% less eutrophicated areas were obtained compared to the areas observed with NDVI and NDAVI. This may be due to the lighting conditions and water brightness, which could overestimate the indices [61,63,64]. Therefore, there is a need to use several vegetation indices for better predictions of vegetation (and eutrophication) patterns and dynamics, as these patterns cannot be determined exclusively by relying on one index due to their different limitations [61].

By checking the average annual CHL and TUR concentrations, it can be seen that they all correspond to the worst trophic state and turbidity levels considered by the LTBAA. These results are consistent with the narrative presented above, which indicates that the lake’s waters are becoming increasingly impacted and that the extent of the contamination is spreading. The analysis of the concentrations showed a decrease in the dry season and an increase in the rainy season, which may be due to nutrient flushing by rain and sediment resuspension, which reiterates the impact that the land cover of the basin has on the water quality of the lake. Once again, we emphasize the importance of quantifying the effects of land cover changes on hydrological variables such as runoff and discharge, as stressed by Srivastava et al. [13] and Wei et al. [12]. In Figure 6, the spatial behavior of TUR appears to be opposite to that of CHL. In areas where the CHL concentration is higher, TUR exhibits low values and vice versa. This could arise because the broad wavelength spectral data on the satellite do not allow proper discrimination of CHL on waters with high suspended sediments, and more so in highly turbid and eutrophic waters [75]. Also, a more significant increase in turbidity concentration than in chlorophyll was identified.

Cohana Bay, located at the outlet of the basin, is known to have waters among the most eutrophic of the lake due to the contamination caused by the upstream urban areas and agricultural activities [71]. This is due to nutrient enrichment and other pollutants from domestic and industrial wastewater from the city of El Alto and nearby urban centers, as well as mining activities scattered throughout the study area [46]. This is also supported by an analysis of water-quality indicators in the rivers of the basin [38] that concluded that the majority of environmental driving forces and pressures are spatially located in the upper parts of the basin, corresponding to mining and urban/industrial areas.

By complementing the trajectory analysis with the analysis of SAVI, NDVI, NDAVI, CHL, and TUR in Lake Titicaca, a better understanding of the extent of existing water-quality problems can be acquired. This driver–pressure–impact relation emphasizes the crucial need to monitor the sources of contamination to control them and improve the water quality not only of the rivers in the basin, but also of the lake. Assessing the impacts of on-site LULC changes, on- and off-site water quality is essential for formulating rational management and development strategies for water resources [12].

Based on this study, an increasing behavior of eutrophication area in the lake was seen at a ratio of approximately 1:4.8 when the urban area increased. An average increase rate of 8 km²/year of urban growth was detected, and a rate of 1.69 km²/year of eutrophicated areas was calculated. This indicates that for the 12-year study period for every 4.8 km² of urban built-up area, the extent of eutrophicated areas in Lake Titicaca’s shores increased by 1 km² on average. Considering that El Alto has such rapid urban growth, if water resources continue to be managed in the same way, Lake Titicaca and its biota do not have a promising outlook. Special emphasis is placed on the fact that Lake Titicaca has high residence times [51] and constitutes the most important water resource of the Bolivian Altiplano [41,46,71], which means that the condition of its waters directly impacts the rural indigenous population and the economy of the country.

Water management is complex due to the interacting subsystems of geological, climatological, ecological, and human nature that are enclosed in it [35]. Thus, finding solutions for water management calls for a deep understanding of its dynamics. Having access to information on the state of the natural resources and what management choices affect and imply, helps with selecting the most suitable choices [35]. These management decisions must be assisted by embedding them in policy and governance on different scales [36]. The outcomes of this research contribute to understanding the impacts of the urban growth driver on water quality, which can assist managers in developing policy in the direction of long-term sustainability.

5. Conclusions

This research aimed to assess the spatial and temporal dynamics of LC change in the KRB, and to explore driver–impact relations on the water quality of Lake Titicaca. A novel approach based on understanding the relation between on-site drivers and their off-site effect through remote-sensing data was followed.

A 123% expansion of urban areas was estimated in 2018 compared to 2006, which corresponded with the increasing trend found in data (p < 0.05) and an increasing rate of 8.023 km²/year for the 12-year study period. This finding evidenced the rapid growth of El Alto city mentioned so often in literature [38,39,41,71]. The Mann–Kendall test perceived trends only in the built-up class, which may be due to the classification deficiency found in Table 3, especially in the vegetation and shrubs/grass classes. In future research, this can be improved by employing the use of in situ data for land cover classification, improving the spatial resolution of satellite images, or by combining multiple remote-sensing data.

Furthermore, by evaluating the behavior of the extent of eutrophication at the basin outlet by analyzing SAVI, NDVI, and NDAVI, it was possible to detect unusually high index values. Increasing trends were found for the vegetation indices based on the Mann–Kendall test with fluctuating magnitudes of change based on Sen’s slopes estimator. A ratio of approximately 1:4.8 among eutrophicated areas downstream versus urbanized areas upstream was calculated. Similarly, when comparing average annual concentrations of CHL and TUR at the lake, an increase of around 40% and 155% was calculated from 2018 to 2016, respectively. These findings, together with findings reported in other studies [38,39,41], showcase the driver–impact relation that the organic waste coming from the urban areas has on the water quality of rivers and the lake. However, the impact of the increasing agricultural areas should also be considered in the future, as this LULC also alters the hydrological and water-quality characteristics [13]. Furthermore, for future research, the impact of climate change on the climatic conditions of the basin should be evaluated [49]. A thorough understanding of all the drivers’ impacts must be achieved so an adjustment in the management of the system in the direction of long-term sustainability can be made [35].

In addition, the basin is located within the jurisdiction of several municipalities and even a binational organization between Perú and Bolivia, which makes it challenging to assign responsibilities and coordinate the management of resources within the region. The institutions that have a presence in the watershed have been allocating resources to end-of-pipe solutions for approximately 15 years. These findings reveal the impact that poor management of solid waste and domestic and industrial wastewater is having on the basin, which is already causing health problems in the population and the aquatic biota.

The trajectory map obtained in this study can be used as a visual tool for the prioritization of areas that are vulnerable to contamination and to better understand the source–impact relation in the river basin. The ministry has shared its interest on the potential of using this visual analysis for explaining, sharing, and consulting with stakeholders future measures to deal with the water challenges in the basin. Going forward, it could be used for the selection of potential water-quality monitoring sites in the rivers, to establish measures for the reduction, mitigation, and treatment of pollution. Eventually, it could be combined with research such as the one developed by Duquesne et al. [43], where a model is used for predicting the impact of human activities on Lake Titicaca’s ecosystem, and therefore, can potentially assist local authorities in managing the lake and the basins surrounding it.

Although this study generated useful information for decision-makers, it is subject to some limitations. Opportunities for improvement were identified in the classification of vegetation classes. When analyzing the error matrix in Table 3, the classes with lower accuracies were vegetation and shrubs/grass. This could be due to close spectral signatures between the classes [57], or because of errors when assigning the classes of the polygons used for the accuracy assessment. Alternatively, a different number of classes may need to be defined to achieve a more accurate classification. To ensure a better classification of the vegetated areas, historical in situ data to calibrate and validate should be used. While Landsat 7 satellite images meet the demand of our study, classification errors could be reduced by combining remote-sensing images with different time and spatial resolutions, such as Sentinel images. Also, a more frequent analysis of LC changes, vegetation indices, and water quality would provide more robust and representative data and trends.

In addition, a higher spatial resolution should be used to properly identify smaller towns and roads, which ought to be included in urban area growth quantification. Similarly, a thorough analysis of the lake’s water quality through in situ measurements and other remote-sensing products should be implemented. Also, future land cover classifications should consider the vegetation (indices) response to precipitation to improve accuracy of analysis.

To ensure the adoption of remote-sensing water-quality monitoring in the long term, it is advisable to make comparisons with in situ data to ensure that the algorithms and band combination used in the retrieval of CHL and TUR are appropriate for this location.

Similarly, a determination of the Optical Water Type present in Lake Titicaca could be useful to get a more comprehensive overview of the optical properties of its waters and take it into account for future water-quality indicators to be retrieved using remote sensing.

Author Contributions

Conceptualization, A.B. and A.A.; Methodology, A.B., A.A., I.R. and A.v.G.; Analysis, A.B. and A.A.; Resources, A.B., A.A. and I.R.; Writing—original draft preparation, A.B. and A.A.; Writing—review and editing, A.B., A.A., I.R. and A.v.G.; Supervision, I.R. and A.v.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research has been supported by the EU H2020 project Water-ForCE (grant no. 101004186) and the VLIR-UOS, IUC 2017 Phase 1 UCB-B, https://www.vliruos.be/en/projects/project/22?pid=3607, accessed on 5 February 2022.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that support the findings of this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Municipalities that encompass the KRB and main rivers.

Figure 2. Land cover maps with a three-year time step. Images corresponding to dry season (a–e) and wet season (f–j). Unclassified patches correspond to masked clouds.

Figure 3. Trajectory map of the Katari River Basin from 2006 to 2018.

Figure 4. NDVI, SAVI, and NDAVI areas in Lake Titicaca for dry and wet periods.

Figure 5. Increase of eutrophication extent in Lake Titicaca measured by NDVI for the study period. Left: dry season, right: wet season.

Figure 6. Distribution of CHL and TUR concentrations.

Table 1

Summary of data downloaded for this study.

Data	Data Description	Data Source
Satellite Images	Landsat 7, level 1 products (Path 1, Rows 71, and 72) Seasonal images with a 3 year time step (study period 2006–2018)	United States Geological Survey last accessed on 28 June 2021https://earthexplorer.usgs.gov/
Digital Elevation Model (DEM)	30 m spatial resolution, obtained from ASTER Global Digital Elevation Model (Tile: S17W069)	GeoBolivia last accessed on 31 May 2021https://geo.gob.bo/portal/?Descargas-Aster-DEM
Chlorophyll-a data (CHL)	Sentinel 2 derived products with a 20 m spatial resolution (Tile: 19KEB) from 2015–2018. CHL data in mg/m3 and TUR data in FNU	Terrascope last accessed on 27 January 2022https://terrascope.be/en
Turbidity data (TUR)

Table 2

Land cover classes in area (km²) and % obtained from LC classification for dry season.

Class	2006		2009		2012		2015		2018
	Area	%	Area	%	Area	%	Area	%	Area	%
Water	22.54	0.7	3.49	0.1	20.52	0.7	11.71	0.4	6.97	0.2
Built-up	79.27	2.6	109.56	3.6	135.29	4.5	150.77	5.0	176.49	5.9
Vegetation	131.50	4.4	129.71	4.3	272.83	9.1	126.64	4.2	77.87	2.6
Barren land	2470.65	82.2	2660.68	88.5	2333.91	77.6	2237.38	74.4	2463.83	82.0
Shrubs/grass	299.07	9.9	99.67	3.3	238.11	7.9	476.98	15.9	263.29	8.8
Snow cover	3.23	0.1	3.15	0.1	5.61	0.2	2.82	0.1	17.79	0.6

Table 3

Error matrix of the land cover classification from 2006 to 2018.

Year	Land Cover	Water	Built-Up	Vegetation	Barren Land	Shrubs/Grass	Snow Cover
2006	PA	100	95.81	76.26	99.82	76.42	100
	UA	100	100	86.03	95.84	96.41	100
	OA	95.62
	Kappa	0.87
2009	PA	100	53.6	55.45	99.12	62.05	100
	UA	100	100	77.95	91.56	97.3	100
	OA	91.49
	Kappa	0.68
2012	PA	100	100	98.54	99.33	35.28	100
	UA	100	100	92.34	82.53	96.3	100
	OA	85.45
	Kappa	0.69
2015	PA	100	97.79	100	97.71	68.08	100
	UA	60.86	100	22.58	95.57	89.25	100
	OA	91.58
	Kappa	0.8
2018	PA	100	100	80.83	99.83	65.22	79.28
	UA	100	100	68.11	94.63	99.03	98.14
	OA	94.67
	Kappa	0.84

Table 4

Mean CHL and TUR concentrations throughout the years (* only one image was available for this year).

Year	2015	2016	2017	2018
CHL (mg/m3)	119.76 *	39.57	49.66	55.45
TUR (FNU)	12.18	17.18	24.87	43.83

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/w14071021/s1, Table S1: Description of satellite images used; Table S2: Mann-Kendall and Sen’s slopes results for land cover classes; Table S3: Extent of high NDVI, NDAVI, and SAVI areas in km² for the study period, and Table S4: Mann-Kendall and Sen’s slopes results for SAVI, NDVI, and NDAVI.

References

1. Mallupattu, P.K.; Sreenivasula Reddy, J.R. Analysis of Land Use/Land Cover Changes Using Remote Sensing Data and GIS at an Urban Area, Tirupati, India. Sci. World J.; 2013; 2013, 268623. [DOI: https://dx.doi.org/10.1155/2013/268623] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23781152]

2. Bonansea, M.; Bazán, R.; Germán, A.; Ferral, A.; Beltramone, G.; Cossavella, A.; Pinotti, L. Assessing Land Use and Land Cover Change in Los Molinos Reservoir Watershed and the Effect on the Reservoir Water Quality. J. South Am. Earth Sci.; 2021; 108, 103243. [DOI: https://dx.doi.org/10.1016/j.jsames.2021.103243]

3. Huang, J.; Zhan, J.; Yan, H.; Wu, F.; Deng, X. Evaluation of the Impacts of Land Use on Water Quality: A Case Study in The Chaohu Lake Basin. Sci. World J.; 2013; 2013, 329187. [DOI: https://dx.doi.org/10.1155/2013/329187] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23970833]

4. Wilson, C. Land Use/Land Cover Water Quality Nexus: Quantifying Anthropogenic Influences on Surface Water Quality. Environ. Monit. Assess.; 2015; 187, 424. [DOI: https://dx.doi.org/10.1007/s10661-015-4666-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26065891]

5. Wilson, C.; Weng, Q. Assessing Surface Water Quality and Its Relation with Urban Land Cover Changes in the Lake Calumet Area, Greater Chicago. Environ. Manag.; 2010; 45, pp. 1096-1111. [DOI: https://dx.doi.org/10.1007/s00267-010-9482-6]

6. Hasmadi, M.; HZ, P.; MF, S. Evaluating Supervised and Unsupervised Techniques for Land Cover Mapping Using Remote Sensing Data. Geogr. Malays. J. Soc. Space; 2009; 5, pp. 1-10.

7. Jeevalakshmi, D.; Reddy, S.N.; Manikiam, B. Land Cover Classification Based on NDVI Using LANDSAT8 Time Series: A Case Study Tirupati Region. Proceedings of the 2016 International Conference on Communication and Signal Processing (ICCSP); Melmaruvathur, India, 6–8 April 2016; pp. 1332-1335.

8. Zaidi, S.M.; Akbari, A.; Samah, A.A.; Kong, N.S.; Gisen, J.I.A. Landsat-5 Time Series Analysis for Land Use/Land Cover Change Detection Using NDVI and Semi-Supervised Classification Techniques. Pol. J. Environ. Stud.; 2017; 26, pp. 2833-2840. [DOI: https://dx.doi.org/10.15244/pjoes/68878]

9. Kc, A.; Wagle, N.; Acharya, T.D. Spatiotemporal Analysis of Land Cover and the Effects on Ecosystem Service Values in Rupandehi, Nepal from 2005 to 2020. ISPRS Int. J. Geo-Inf.; 2021; 10, 635. [DOI: https://dx.doi.org/10.3390/ijgi10100635]

10. Gong, J.; Sui, H.; Ma, G.; Zhou, Q. A Review of Multi-Temporal Remote Sensing Data Change Detection Algorithms. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2008; 37, pp. 757-762.

11. Kuma, H.G.; Feyessa, F.F.; Demissie, T.A. Land-Use/Land-Cover Changes and Implications in Southern Ethiopia: Evidence from Remote Sensing and Informants. Heliyon; 2022; 8, e09071. [DOI: https://dx.doi.org/10.1016/j.heliyon.2022.e09071]

12. Wei, P.; Chen, S.; Wu, M.; Deng, Y.; Xu, H.; Jia, Y.; Liu, F. Using the InVEST Model to Assess the Impacts of Climate and Land Use Changes on Water Yield in the Upstream Regions of the Shule River Basin. Water; 2021; 13, 1250. [DOI: https://dx.doi.org/10.3390/w13091250]

13. Srivastava, A.; Kumari, N.; Maza, M. Hydrological Response to Agricultural Land Use Heterogeneity Using Variable Infiltration Capacity Model. Water Resour Manag.; 2020; 34, pp. 3779-3794. [DOI: https://dx.doi.org/10.1007/s11269-020-02630-4]

14. Garg, V.; Nikam, B.R.; Thakur, P.K.; Aggarwal, S.P.; Gupta, P.K.; Srivastav, S.K. Human-Induced Land Use Land Cover Change and Its Impact on Hydrology. HydroResearch; 2019; 1, pp. 48-56. [DOI: https://dx.doi.org/10.1016/j.hydres.2019.06.001]

15. Getachew, H.; Melesse, A. The Impact of Land Use Change on the Hydrology of the Angereb Watershed, Ethiopia. Int. J. Water Sci.; 2012; 1, 4. [DOI: https://dx.doi.org/10.5772/56266]

16. Koneti, S.; Sunkara, S.L.; Roy, P.S. Hydrological Modeling with Respect to Impact of Land-Use and Land-Cover Change on the Runoff Dynamics in Godavari River Basin Using the HEC-HMS Model. ISPRS Int. J. Geo-Inf.; 2018; 7, 206. [DOI: https://dx.doi.org/10.3390/ijgi7060206]

17. Alvarenga, L.A.; de Mello, C.R.; Colombo, A.; Cuartas, L.A.; Bowling, L.C. Assessment of Land Cover Change on the Hydrology of a Brazilian Headwater Watershed Using the Distributed Hydrology-Soil-Vegetation Model. Catena; 2016; 143, pp. 7-17. [DOI: https://dx.doi.org/10.1016/j.catena.2016.04.001]

18. Boori, M.; Voženílek, V. Remote Sensing and Land Use/Land Cover Trajectories. J. Geophys. Remote Sens.; 2014; 3, 107. [DOI: https://dx.doi.org/10.4172/2169-0049.1000123]

19. Zioti, F.; Ferreira, K.R.; Queiroz, G.R.; Neves, A.K.; Carlos, F.M.; Souza, F.C.; Santos, L.A.; Simoes, R.E.O. A Platform for Land Use and Land Cover Data Integration and Trajectory Analysis. Int. J. Appl. Earth Obs. Geoinf.; 2022; 106, 102655. [DOI: https://dx.doi.org/10.1016/j.jag.2021.102655]

20. De Keukelaere, L.; Knaeps, E. Terrascope Sentinel-2 Algorithm Theoretical Base Document S2—WATER QUALITY—V100. 2021; Available online: https://docs.terrascope.be/DataProducts/Sentinel-2/references/VITO_S2_ATBD_S2_WATER_QUALITY_V100.pdf (accessed on 27 December 2021).

21. Ozbay, G.; Fan, C.; Yang, Z. Relationship between Land Use and Water Quality and Its Assessment Using Hyperspectral Remote Sensing in Mid- Atlantic Estuaries. Water Quality; Intech: London, UK, 2017; 428. ISBN 978-953-51-2881-6

22. Ouma, Y.O.; Noor, K.; Herbert, K. Modelling Reservoir Chlorophyll-a, TSS, and Turbidity Using Sentinel-2A MSI and Landsat-8 OLI Satellite Sensors with Empirical Multivariate Regression. J. Sens.; 2020; 2020, e8858408. [DOI: https://dx.doi.org/10.1155/2020/8858408]

23. Caballero, I.; Stumpf, R.P.; Meredith, A. Preliminary Assessment of Turbidity and Chlorophyll Impact on Bathymetry Derived from Sentinel-2A and Sentinel-3A Satellites in South Florida. Remote Sens.; 2019; 11, 645. [DOI: https://dx.doi.org/10.3390/rs11060645]

24. Katlane, R.; Dupouy, C.; Kilani, B.E.; Berges, J.C. Estimation of Chlorophyll and Turbidity Using Sentinel 2A and EO1 Data in Kneiss Archipelago Gulf of Gabes, Tunisia. Int. J. Geosci.; 2020; 11, 708. [DOI: https://dx.doi.org/10.4236/ijg.2020.1110035]

25. Kuhn, C.; de Matos Valerio, A.; Ward, N.; Loken, L.; Sawakuchi, H.O.; Kampel, M.; Richey, J.; Stadler, P.; Crawford, J.; Striegl, R. et al. Performance of Landsat-8 and Sentinel-2 Surface Reflectance Products for River Remote Sensing Retrievals of Chlorophyll-a and Turbidity. Remote Sens. Environ.; 2019; 224, pp. 104-118. [DOI: https://dx.doi.org/10.1016/j.rse.2019.01.023]

26. Watanabe, F.; Alcântara, E.; Rodrigues, T.; Rotta, L.; Bernardo, N.; Imai, N. Remote Sensing of the Chlorophyll-a Based on OLI/Landsat-8 and MSI/Sentinel-2A (Barra Bonita Reservoir, Brazil). An. Acad. Bras. Ciênc.; 2017; 90, pp. 1987-2000. [DOI: https://dx.doi.org/10.1590/0001-3765201720170125] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28876398]

27. Pahlevan, N.; Smith, B.; Schalles, J.; Binding, C.; Cao, Z.; Ma, R.; Alikas, K.; Kangro, K.; Gurlin, D.; Hà, N. et al. Seamless Retrievals of Chlorophyll-a from Sentinel-2 (MSI) and Sentinel-3 (OLCI) in Inland and Coastal Waters: A Machine-Learning Approach. Remote Sens. Environ.; 2020; 240, 111604. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111604]

28. Garg, V.; Aggarwal, S.P.; Chauhan, P. Changes in Turbidity along Ganga River Using Sentinel-2 Satellite Data during Lockdown Associated with COVID-19. Geomat. Nat. Hazards Risk; 2020; 11, pp. 1175-1195. [DOI: https://dx.doi.org/10.1080/19475705.2020.1782482]

29. Ansper, A.; Alikas, K. Retrieval of Chlorophyll a from Sentinel-2 MSI Data for the European Union Water Framework Directive Reporting Purposes. Remote Sens.; 2019; 11, 64. [DOI: https://dx.doi.org/10.3390/rs11010064]

30. Patel, H.; Vashi, R.T. Chapter 6—Use of Naturally Prepared Coagulants for the Treatment of Wastewater from Dyeing Mills. Characterization and Treatment of Textile Wastewater; Patel, H.; Vashi, R.T. Elsevier: Boston, MA, USA, 2015; pp. 147-158. ISBN 978-0-12-802326-6

31. Jaskuła, J.; Sojka, M. Assessing Spectral Indices for Detecting Vegetative Overgrowth of Reservoirs. Pol. J. Environ. Stud.; 2019; 28, pp. 4199-4211. [DOI: https://dx.doi.org/10.15244/pjoes/98994]

32. Villa, P.; Laini, A.; Bresciani, M.; Bolpagni, R. A Remote Sensing Approach to Monitor the Conservation Status of Lacustrine Phragmites Australis Beds. Wetl. Ecol. Manag.; 2013; 21, pp. 399-416. [DOI: https://dx.doi.org/10.1007/s11273-013-9311-9]

33. Villa, P.; Bresciani, M.; Braga, F.; Bolpagni, R. Comparative Assessment of Broadband Vegetation Indices Over Aquatic Vegetation. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2014; 7, pp. 3117-3127. [DOI: https://dx.doi.org/10.1109/JSTARS.2014.2315718]

34. Miheretu, B.A.; Yimer, A.A. Land Use/Land Cover Changes and Their Environmental Implications in the Gelana Sub-Watershed of Northern Highlands of Ethiopia. Env. Syst Res; 2017; 6, 7. [DOI: https://dx.doi.org/10.1186/s40068-017-0084-7]

35. Keesstra, S.; Mol, G.; De Leeuw, J.; Okx, J.; Molenaar, C.; De Cleen, M.; Visser, S. Soil-Related Sustainable Development Goals: Four Concepts to Make Land Degradation Neutrality and Restoration Work. Land; 2018; 7, 133. [DOI: https://dx.doi.org/10.3390/land7040133]

36. Keesstra, S.; Sannigrahi, S.; López-Vicente, M.; Pulido, M.; Novara, A.; Visser, S.; Kalantari, Z. The Role of Soils in Regulation and Provision of Blue and Green Water. Philos. Trans. R. Soc. B Biol. Sci.; 2021; 376, 20200175. [DOI: https://dx.doi.org/10.1098/rstb.2020.0175] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34365828]

37. Gorgoglione, A.; Gregorio, J.; Ríos, A.; Alonso, J.; Chreties, C.; Fossati, M. Influence of Land Use/Land Cover on Surface-Water Quality of Santa Lucía River, Uruguay. Sustainability; 2020; 12, 4692. [DOI: https://dx.doi.org/10.3390/su12114692]

38. Agramont, A.; van Cauwenbergh, N.; van Griesven, A.; Craps, M. Integrating Spatial and Social Characteristics in the DPSIR Framework for the Sustainable Management of River Basins: Case Study of the Katari River Basin, Bolivia. Water Int.; 2021; 47, pp. 8-29. [DOI: https://dx.doi.org/10.1080/02508060.2021.1997021]

39. Chudnoff, S. A Water Quality Assessment of the Rio Katari River and Its Principle Tributaries, Bolivia. 2009; Available online: https://digitalrepository.unm.edu/wr_sp/127/ (accessed on 5 February 2022).

40. CGEMinisterio de Energía y MetalurgiaMMAyAMinisterio de Relaciones ExterioresGobierno Autónomo Departamental de la PazGobiernos Autónomos Municipales de El Alto, Viacha, Laja, Pucarani y Puerto PérezEPSASEMALT. FPS Informe de Auditoría Sobre El Desempeño Ambiental Respecto de La Contaminación Hídrica En La Cuenca Del Río Katari y La Bahía de Cohana.Contraloría General del Estado (CGE), Bolivia. 2014; Available online: https://www.contraloria.gob.bo/portal/Auditor%C3%ADa/Auditor%C3%ADaAmbiental.aspx (accessed on 11 November 2021).

41. Archundia, D.; Duwig, C.; Spadini, L.; Uzu, G.; Guédron, S.; Morel, M.C.; Cortez, R.; Ramos Ramos, O.; Chincheros, J.; Martins, J.M.F. How Uncontrolled Urban Expansion Increases the Contamination of the Titicaca Lake Basin (El Alto, La Paz, Bolivia). Water Air Soil Pollut.; 2016; 228, 44. [DOI: https://dx.doi.org/10.1007/s11270-016-3217-0]

42. Pahl-Wostl, C. A Methodological Framework for Empirical Analysis. Water Governance in the Face of Global Change: From Understanding to Transformation; Pahl-Wostl, C. Water Governance—Concepts, Methods, and Practice Springer International Publishing: Cham, Switzerland, 2015; pp. 181-201. ISBN 978-3-319-21855-7

43. Duquesne, F.; Vallaeys, V.; Vidaurre, P.J.; Hanert, E. A Coupled Ecohydrodynamic Model to Predict Algal Blooms in Lake Titicaca. Ecol. Model.; 2021; 440, 109418. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2020.109418]

44. Agramont, A.; Craps, M.; Balderrama, M.; Huysmans, M. Transdisciplinary Learning Communities to Involve Vulnerable Social Groups in Solving Complex Water-Related Problems in Bolivia. Water; 2019; 11, 385. [DOI: https://dx.doi.org/10.3390/w11020385]

45. Jhonatan, E.; Vallejos, A.; Saavedra, O.; Escalera, A. Evaluación de la precipitación distribuida en la cuenca katari basado en tecnología satelital y productos derivados. Investig. Desarro.; 2018; 18, pp. 35-51. [DOI: https://dx.doi.org/10.23881/idupbo.018.1-3i]

46. Molina, C.; Lazzaro, X.; Guedron, S.; Acha, C.D. Contaminación de La Bahía de Cohana, Lago Titicaca (Bolivia): Desafíos y Oportunidades Para Promover Su Recuperación Pollution at Cohana Bay, Lake Titicaca (Bolivia): Challenges and Opportunities to Promote Its Recovery. Ecol. Boliv.; 2017; 52, pp. 65-76.

47. BID Programa de Saneamiento Del Lago Titicaca (Cuenca Katari, Bahía Cohana). 2016; Available online: https://ewsdata.rightsindevelopment.org/files/documents/18/IADB-BO-L1118_PbsvUa3.pdf (accessed on 11 November 2021).

48. Arbona, J.M.; Kohl, B. La Paz–El Alto. Cities; 2004; 21, pp. 255-265. [DOI: https://dx.doi.org/10.1016/j.cities.2004.02.004]

49. Buxton, N.; Escobar, M.; Purkey, D.; Lima, N. Water Scarcity, Climate Change and Bolivia: Planning for Climate Uncertainties; SEI Discussion Brief Stockholm Environment Institute: Stockholm, Sweden, 2013.

50. MMAyA Plan Director de la Cuenca Katari. 2010; Available online: https://www.academia.edu/27074895/Plan_director_katari (accessed on 28 June 2021).

51. Rieckermann, J.; Daebel, H.; Ronteltap, M.; Bernauer, T. Assessing the Performance of International Water Management at Lake Titicaca. Aquat. Sci.; 2006; 68, pp. 502-516. [DOI: https://dx.doi.org/10.1007/s00027-006-0863-0]

52. GISGeography. What Is Atmospheric Correction in Remote Sensing?. Available online: https://gisgeography.com/atmospheric-correction/ (accessed on 17 March 2022).

53. Congedo, L. Semi-Automatic Classification Plugin: A Python Tool for the Download and Processing of Remote Sensing Images in QGIS. J. Open Source Softw.; 2021; 6, 3172. [DOI: https://dx.doi.org/10.21105/joss.03172]

54. Li, X.C.; Liu, L.L.; Huang, L.K. Comparison of several remote sensing image classification methods based on envi. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2020; 42, pp. 605-611. [DOI: https://dx.doi.org/10.5194/isprs-archives-XLII-3-W10-605-2020]

55. Abbas, A.; Minallh, N.; Ahmad, N.; Abdur Rehman Abid, S.; Khan, M. K-Means and ISODATA Clustering Algorithms for Landcover Classification Using Remote Sensing. Sindh Univ. Res. J.; 2016; 48, pp. 315-318.

56. Vimala, R.; Marimuthu, A.; Venkateswaran, S.; Poongodi, R. Unsupervised ISODATA Algorithm Classification Used in the Landsat Image for Predicting the Expansion of Salem Urban, Tamil Nadu. Indian J. Sci. Technol.; 2020; 13, pp. 1619-1629. [DOI: https://dx.doi.org/10.17485/IJST/v13i16.271]

57. Manandhar, R.; Odeh, I.O.A.; Ancev, T. Improving the Accuracy of Land Use and Land Cover Classification of Landsat Data Using Post-Classification Enhancement. Remote Sens.; 2009; 1, pp. 330-344. [DOI: https://dx.doi.org/10.3390/rs1030330]

58. USGS. Landsat Normalized Difference Vegetation Index. Available online: https://www.usgs.gov/core-science-systems/nli/landsat/landsat-normalized-difference-vegetation-index?qt-science_support_page_related_con=0#qt-science_support_page_related_con (accessed on 19 June 2021).

59. Zhou, Q.; Li, B.; Kurban, A. Trajectory Analysis of Land Cover Change in Arid Environment of China. Int. J. Remote Sens.; 2008; 29, pp. 1093-1107. [DOI: https://dx.doi.org/10.1080/01431160701355256]

60. Hussain, M.M.; Mahmud, I. PyMannKendall: A Python Package for Non Parametric Mann Kendall Family of Trend Tests. J. Open Source Softw.; 2019; 4, 1556. [DOI: https://dx.doi.org/10.21105/joss.01556]

61. Kumari, N.; Srivastava, A.; Dumka, U.C. A Long-Term Spatiotemporal Analysis of Vegetation Greenness over the Himalayan Region Using Google Earth Engine. Climate; 2021; 9, 109. [DOI: https://dx.doi.org/10.3390/cli9070109]

62. Xiao, Q.; Zhang, M.; Hu, Z.; Gao, Y.; Hu, C.; Liu, C.; Liu, S.; Zhang, Z.; Zhao, J.; Xiao, W. et al. Spatial Variations of Methane Emission in a Large Shallow Eutrophic Lake in Subtropical Climate. J. Geophys. Res. Biogeosci.; 2017; 122, pp. 1597-1614. [DOI: https://dx.doi.org/10.1002/2017JG003805]

63. Martín-Ortega, P.; García-Montero, L.G.; Sibelet, N. Temporal Patterns in Illumination Conditions and Its Effect on Vegetation Indices Using Landsat on Google Earth Engine. Remote Sens.; 2020; 12, 211. [DOI: https://dx.doi.org/10.3390/rs12020211]

64. Xue, J.; Su, B. Significant Remote Sensing Vegetation Indices: A Review of Developments and Applications. J. Sens.; 2017; 2017, e1353691. [DOI: https://dx.doi.org/10.1155/2017/1353691]

65. USGS. Landsat Soil Adjusted Vegetation Index. U.S. Geological Survey. Available online: https://www.usgs.gov/landsat-missions/landsat-soil-adjusted-vegetation-index (accessed on 15 March 2022).

66. Villa, P.; Mousivand, A.; Bresciani, M. Aquatic Vegetation Indices Assessment through Radiative Transfer Modeling and Linear Mixture Simulation. Int. J. Appl. Earth Obs. Geoinf.; 2014; 30, pp. 113-127. [DOI: https://dx.doi.org/10.1016/j.jag.2014.01.017]

67. Anderson, J.R.; Hardy, E.E.; Roach, J.T.; Witmer, R.E. A Land Use and Land Cover Classification System for Use with Remote Sensor Data; Professional Paper US Government Printing Office: Washington, DC, USA, 1976.

68. LTBAA Parámetros Indicadores y Criterios de Evaluación de La Situación de Las Aguas Del Lago. Available online: http://www.alt-perubolivia.org/web_lago/WEB_LT/Finales/accion/accion_2.htm#213 (accessed on 28 December 2021).

69. Zoungrana, B.J.-B.; Conrad, C.; Amekudzi, L.K.; Thiel, M.; Da, E.D. Land Use/Cover Response to Rainfall Variability: A Comparing Analysis between NDVI and EVI in the Southwest of Burkina Faso. Climate; 2015; 3, pp. 63-77. [DOI: https://dx.doi.org/10.3390/cli3010063]

70. Casas, C.; Escobar, J. Diseño de redes de monitoreo apoyadas por herramientas sig y modelación geoespacial. Aqua-LAC; 2011; 3, pp. 18-25. [DOI: https://dx.doi.org/10.29104/phi-aqualac/2011-v3-1-03]

71. Duwig, C.; Archundia, D.; Lehembre, F.; Spadini, L.; Morel, M.C.; Uzu, G.; Chincheros, J.; Cortez, R.; Martins, J.M.F. Impacts of Anthropogenic Activities on the Contamination of a Sub Watershed of Lake Titicaca. Are Antibiotics a Concern in the Bolivian Altiplano?. Procedia Earth Planet. Sci.; 2014; 10, pp. 370-375. [DOI: https://dx.doi.org/10.1016/j.proeps.2014.08.062]

72. Saravanan, S.; Jegankumar, R.; Selvaraj, A.; Jacinth Jennifer, J.; Parthasarathy, K.S.S. Chapter 20—Utility of Landsat Data for Assessing Mangrove Degradation in Muthupet Lagoon, South India. Coastal Zone Management; Ramkumar, M.; James, R.A.; Menier, D.; Kumaraswamy, K. Elsevier: Amsterdam, The Netherlands, 2019; pp. 471-484. ISBN 978-0-12-814350-6

73. Viana, C.M.; Oliveira, S.; Oliveira, S.C.; Rocha, J. Land Use/Land Cover Change Detection and Urban Sprawl Analysis. Spatial Modeling in GIS and R for Earth and Environmental Sciences; Pourghasemi, H.R.; Gokceoglu, C. Elsevier: Amsterdam, The Netherlands, 2019; pp. 621-651. ISBN 978-0-12-815226-3

74. Vargas-Cuentas, N.I.; Roman-Gonzalez, A. Analysis of Harmful Algal Blooms in Lake Titicaca Using Remote Sensing. 21 October 2019. Available online: https://www.researchgate.net/publication/337945152_Analysis_of_harmful_algal_blooms_in_Lake_Titicaca_using_remote_sensing (accessed on 14 December 2021).

75. Gholizadeh, M.H.; Melesse, A.M.; Reddi, L. A Comprehensive Review on Water Quality Parameters Estimation Using Remote Sensing Techniques. Sensors; 2016; 16, 1298. [DOI: https://dx.doi.org/10.3390/s16081298]