Mapping hydropower expansion and cash crop dynamics in Colombia using Landsat time series

Abstract

Satellite-based monitoring provides insights on synergies and trade-offs between energy production and land use and land cover (LULC) changes (LULCC) around hydropower dams. Using Landsat data, we mapped LULCC related to hydropower expansion to understand change dynamics of cash crop production in the Magdalena River basin (Colombian Andes). We leveraged secondary map products and an active learning routine to produce thematically detailed LULC maps, for 2009, 2015, and 2020, including key cash crops (rice and oil palm). Our area-adjusted accuracy assessment revealed high overall and class-specific user’s and producer’s accuracies, exceeding 80% for oil palm, rice, non-agricultural vegetation, wetlands, and grasslands, and accuracies of ∼60% for other temporary and permanent crops. We identified substantial losses of rice (10,6%) and grasslands (25,9%) in El Quimbo due to reservoir flooding. For Hidrosogamoso, there was a loss of wetlands downstream. Our products provide the basis for balancing social, economic, and environmental trade-offs related to conflicting land uses.

Keywords:

• active learning
• Colombian Andes
• Google Earth Engine
• hydropower dams
• land use and land cover (LULC)
• water-energy-food nexus

1. Introduction

Land is a limited resource, underlying increasing pressure for multiple competing uses, including agriculture, extractive activities, energy production, built-up environments, but also for protection of the natural environment, or biodiversity conservation (Munroe and Müller 2020; Xie et al. 2020; Turner et al. 2021; Deteix et al. 2023). Global energy demand has increased rapidly in the last decade (IEA 2022). At the same time, the climate crisis and emerging tipping point discussions (Dakos et al. 2019) have triggered a resurgence of large hydropower projects for renewable energy production, with hotspots of future dam construction clustering in the Global South (Zarfl et al. 2015). Dams induce land system changes, from their construction to the operational phase, altering the access to and the use of land and water resources (Rufin et al. 2019). These changes play different roles downstream and upstream of dams, with drastic changes of local livelihoods (Castro-Diaz et al. 2023) and major ecosystem changes, such as depletion of riparian vegetation and wetlands (Rosenberg et al. 2000). The land surrounding reservoirs is often subject to increased deforestation or agricultural intensification, triggered by irrigation (Rufin et al. 2018; Lohani et al. 2020; Velastegui-Montoya et al. 2020) while at the same time traditional fishing and agriculture are negatively affected (Richter et al. 2010). Already prevailing inequalities may increase, such as disputes over land ownership, water distribution and energy resources (Boelens et al. 2016). Understanding land use and land cover change (LULCC) around dams is hence one of the crucial impacts for monitoring dam-related changes. Such information is also core for disentangling diverse impacts on agricultural land use and productivity, threats to biodiversity and habitats, and carbon sequestration potential (Domínguez and Rivera 2010; Woldemichael et al. 2012; De Oliveira Serrão et al. 2020; De Andrade et al. 2021; Zhang et al. 2022).

The Colombian Andes constitute a prime example of growing pressure on land due to conflicting land uses. The region is a hot-spot of dam construction, while at the same time being a key production area for cash crops. Colombia is planning to increase its current electricity production by 135% until 2050 (UPME 2020) and has a remaining hydropower potential to be explored (Xu et al. 2023), making it a hotspot of future hydropower expansion. Most hydropower dams are concentrated in the Magdalena River basin and the Andes region (Craven et al. 2017; Angarita et al. 2018), where an expansion of the hydropower sector is likely due to the existing infrastructure and the suitable topography of the region (Angarita et al. 2018). The Magdalena basin is home to 70% of the Colombian population (DANE 2022), responsible for 80% of the national energy production, and 70% of national agricultural production (Craven et al. 2017), mostly linked to the production of ‘cash crops’, or ‘export crops’ such as coffee, rice, and oil palm (Agronet 2021). The basin is a biodiversity hotspot (Myers et al. 2000), consisting of wetlands and paramos (Ricaurte et al. 2017; Patiño et al. 2021). Between 2005 and 2010, the Magdalena River basin experienced a loss of 40% of its natural land cover, while agriculture and pasture increased 65% (Restrepo and Escobar 2018). Despite a regrowth of forests in marginal agricultural areas (Rubiano et al. 2017), deforestation linked to the production of cash crops such as coffee and oil palm is one of the major drivers of change in the tropical forests of Colombia (Pendrill et al. 2019b, 2019a; Treanor and Saunders 2021). The growing demand for renewable energy and cash crops (Raihan 2023) requires an improved understanding of the trade-offs between energy production, agricultural production, and the natural ecosystems in this dynamic region.

Land use and land cover (LULC) maps derived from Earth observation data are an important tool to monitor changes and provide timely and standardized information for policies and environmental management (Chowdhury et al. 2021). Data acquired from the Landsat sensor family offers unique data continuity since the 1970s (Chowdhury et al. 2017; Wulder et al. 2019). The Landsat Global Archive Consolidation effort, along with standardized pre-processing schemes and improved processing capabilities now allows using all available data for long and dense time series analyses (Zhu and Woodcock 2014; Wulder et al. 2022). In some world regions, such as in parts of Africa and Latin America, Landsat data availability can be limited, due to persistent cloud cover and limited observation capacity in previous decades (Wulder et al. 2016). As such, integrating time series information across multiple years can be a necessity to obtain gap-free coverage for historic time slices and for separating spectrally similar LULC classes (Griffiths et al. 2013; Müller et al. 2015; Potapov et al. 2021). Multiple studies have investigated the use of aggregating time series information using statistical aggregation functions on distribution and variance - also referred to as spectral temporal metrics (STMs) - in LULC mapping. STM have been successfully used for applications targeting the mapping of agricultural land use across large regions and long time frames (Waldner et al. 2015; Rufin et al. 2022b), involving data-scarce periods and regions (Pflugmacher et al. 2019; Tang et al. 2021; Nill et al. 2022), while maintaining the temporal information that allows for creating thematically detailed maps, e.g. on particular tree species, crop types or management practices (Deines et al. 2017; Hemmerling et al. 2021; Ibrahim et al. 2021).

Multiple map products exist for Colombia, but limitations in thematic and spatial detail, and temporal coverage hamper assessments of LULCC related to hydropower expansion and cash crops. Global maps can provide information on forests (Hansen et al. 2013), water surfaces (Pickens et al. 2020), croplands (Potapov et al. 2021), or specific crop types, such as oil palm (Descals et al. 2021). Such maps can provide a valuable resource for identifying spatial patterns and distribution at the national to global level, but they do not provide the thematic detail for capturing regionally relevant change trajectories and may suffer region-specific inaccuracies (Tulbure et al. 2021). At the national level, the Instituto de Hidrología, Meterología y Estudios Ambientales (IDEAM) provides LULC maps based on the CORINE methodology (IDEAM 2010a). These maps were generated for specific years (2002, 2009, 2012, and 2018) with a minimum mapping unit of 25 ha (IDEAM 2023). Maps created by the SINCHI Institute (SINCHI 2022) and MapBiomas (Mapbiomas 2021) are only available for the Colombian Amazon. González-González et al. (2022) provide a LULC map for the Colombian Andes and Amazon at high spatial resolution for 2018, but provide limited thematic detail in the agriculture domain, which does not allow for in-depth assessments of land change trajectories relevant in the context of dam construction and the dynamics of cash crops in the Magdalena River basin.

In this study, we mapped LULC around two of the major hydropower reservoirs located in the Magdalena River basin (Colombia) with a specific focus on understanding the dynamics of cash crops and dam construction. The Magdalena River basin is particularly cloud-prone due to advective cloud cover from the Andes mountains, resulting in data-scarce periods and regions, which further challenge assessments in these heterogeneous and dynamic landscapes (Ennen et al. 2021). We, therefore, used STMs to aggregate multi-year data into three periods for change assessment, relating to the construction phase of the dams, the operation onset, and the current state after 5 years of operation. The dam licensing inadequately addresses the potential impacts of the dams, considering only adjacent municipalities as the impacted ‘area of influence’ (AOI), using a single map representing the pre-construction phase at coarse spatial resolution and low thematic detail, and does not monitor LULCC for better alignment with other territorial planning (ENEL 2008; ISAGEN 2018). Recognizing that dams can have much broader impacts (Rufin et al. 2019), we select a larger area of interest that allows us to identify potential LULCC at larger spatial scales. We seek to answer three main research questions:

1. How well can we map complex LULC classes in the Magdalena River basin based on sparse Landsat time series?

2. What are the key LULCC trajectories in the reservoir area, in the dam’s area of influence (AOI), and beyond?

3. What are the dynamics of cash crops and how did the dam construction affect these?

2. Data & methods

2.1. Study area

We focus on two hydroelectric dams, El Quimbo, and Hidrosogamoso, in the Magdalena River basin in the Colombian Andes (Figure 1). The region is characterized by a complex topography with gradients from southwest to northeast descending from 5,337 m.a.s.l. in the Colombian Andes to sea level in floodplain areas (Farr et al. 2007). Annual mean temperature and precipitation depend largely on the complex topography, varying from less than 8 to more than 28 °C (IDEAM 2014), and 1000 mm/year in the eastern to 5000 mm/year in the western part of the basin, respectively (IDEAM 2010b). The northern and central parts of the Magdalena basin experience a bimodal annual precipitation cycle with two wet periods (March – May, and October – November) (Restrepo and Escobar 2018). The basin is frequently impacted by El Niño - Southern Oscillation (ENSO) and La Nina events (Hoyos et al. 2013; Poveda et al. 2020). Unique natural vegetation types in the region include the paramos (Patiño et al. 2021) and numerous wetlands (Patino and Estupinan-Suarez 2016; Ricaurte et al. 2017).

Figure 1. Region of interest, (a) El Quimbo dam and (b) Hidrosogamoso, with AOI, Sub-basins upstream and downstream Magdalena River basin, and Andes region.

The agricultural land uses of the Magdalena basin are mainly extensive pastures of low productivity (Etter et al. 2006; Zuluaga et al. 2021), palm oil, coffee, and rice producing (IDEAM 2023). El Quimbo and Hidrosogamoso, located in Huila and Santander departaments, are dominated by different cash crops. In the El Quimbo region, coffee and rice are the most dominant. While rice is mainly produced for the domestic market, coffee targets the domestic and export markets. Hidrosogamoso features extensive oil palm plantations for the international market (Agronet 2021; USDA 2023a).

The El Quimbo and Hidrosogamoso dams were built after the introduction of a new legal framework in 2000 that introduced new guidelines for environmental and social sustainability in dam construction in Colombia. The construction of both dams started in 2009, Hidrosogamoso was commissioned in 2014 and El Quimbo in 2015, and both dams were in full operation in 2015. El Quimbo and Hidrosogamoso are together responsible for 16% of the national energy supply. The license documents of the dams define an AI comprising the municipalities where the dam is installed and the neighboring municipalities (ENEL 2008; ISAGEN 2018). Because downstream impacts of hydropower dams are often neglected due to the frequently long distances between dams and impacted areas, and the less obvious nature of downstream impacts (Baird et al. 2021), we here consider the catchment area that drains towards the dam (dams upstream) and the drainage basins downstream of the dam, including all eastern Andean municipalities in the area downstream (IDEAM 2013) summing up to 3.164.598,82 ha for El Quimbo, and 4.996.476,97 ha for Hidrosogamoso (Figure 1).

2.2. Workflow

The methodology consists of seven main steps (Figure 2), including pre-processing of the Landsat archive into STM, combining secondary map products and an active learning scheme to compile training data for model training, LULC classification, and area-adjusted accuracy assessment.

Figure 2. Workflow of the study.

2.3. Landsat data preprocessing

We used all Landsat 5 TM, 7 ETM+, and 8 OLI Level 2, Collection 2, Tier 1 images available in Google Earth Engine (GEE) between 2007 and 2021. We utilized the Landsat quality assessment (QA) band to mask clouds, cloud shadows and dilated clouds based on the Fmask cloud detection (Qiu et al. 2019). Due to the high cloud cover and limited image availability in parts of the region, we aggregated observations across several years to reach a sufficient number of clear-observation counts for each period of interest. As the focus of our study is on the changes brought about by construction and operation of the hydropower reservoirs, we defined three periods aligning with major developments of these infrastructures, representing before, after and during construction. We integrated the years 2007–2011 for the target year 2009, the years 2013–2017 for the target year 2015 and the years 2018–2022 for the target year 2020. This resulted in 974 images for 2009, 2,075 for 2015, and 1,743 for 2020, leading to varying pixel-level clear observation counts across the study period and region (Figure 3).

Figure 3. Pixel-wise count of clear observations for each aggregation period for (A) Hidrogamoso and (B) El Quimbo.

For each period, we calculated eight spectral-temporal metrics (maximum, minimum, median, standard deviation, 10th, 25th, 75th, and 90th percentile) from eight spectral indices (Table 1), three linear band transformations, and six Landsat spectral bands. We selected these indices as previous studies showed that they have the potential to pick up different aspects of data distributions in the feature space regarding our classes of interest. To incorporate pre-existing knowledge of crop distribution depending on topography (Zuluaga et al. 2021), we included elevation, slope, and aspect from the shuttle radar topography mission (SRTM) data (Farr et al. 2007). Overall, this resulted in a total of 139 features for each period.

Table 1. Spectral indices and authors.

Spectral Index	Application	Reference
Normalized difference vegetation index (NDVI)	Vegetation assessment. Including agricultural vegetation.	Huang et al. (2021); Tucker (1979)
Enhanced vegetation index (EVI)	Vegetation assessment. Reduces saturation effects in canopies with high and low leaf area index (LAI) values, reduces sensitivity to aerosol loads.	Huete et al. (2002)
Soil-adjusted vegetation index (SAVI)	Vegetation assessment. Reduces soil impacts on the vegetation signal.	Huete (1988)
SWIR Ratio	Multispectral surrogate for the Cellulose Absorption Index (CAI) can improve the delineation of dead biomass and soils.	Kowalski et al. (2022)
Normalized difference moisture index (NDMI)	Water surface and wetlands. Sensitive to changes in liquid water in vegetation canopies.	Gao (1996)
Normalized difference burn ratio (NBR)	Land management, applied to grasslands and pastures. Also sensible to water surface. Depicts the relation between NIR and SWIR, highlighting burnt areas.	Jakimow et al. (2018); Key and Benson (2006)
Chlorophyll and carotenoid index (CCI)	Vegetation assessment. Sensitive to phenological changes in evergreen foliage.	Gamon et al. (2016)
Atmospherically resistant vegetation index (ARVI)	Vegetation assessment. Incorporates BLUE to account for high aerosol loads and other atmospheric effects.	Kaufman and Tanré (1992)
Tasseled cap greenness (TCG), brightness (TCB), and wetness (TCW)	Vegetation assessment, incorporates information from the whole feature space including all bands.	Crist (1985); Huang et al. (2002)

See supplementary material (Fig. S1) for the variable importance of the final model.

2.4. Model training and classification

Following the regional characteristics of the land system and the expected change trajectories, we concentrated our analysis on nine classes: rice, oil palm plantations, other permanent crops, other temporary crops, grassland, non-agricultural vegetation, wetlands, water surfaces, and others (Table 2 and Figure 4).

Figure 4. Field photographs with examples of (a) rice cultivation, (b) oil palm cultivation, (c) grasslands, (d) non-agricultural vegetation, (e) water surface, (f) others, (g) other temporary crops, (h) other permanent crops and (i) wetlands were obtained by the authors during fieldwork in 2023.

Table 2. Class definitions.

Class	Definition
Rice (Figure 1c)	Areas with rice cultivation, or in preparation or harvest stage.
Oil palm plantations (Figure 1d)	Vital and actively used palm oil plantations
Grasslands (Figure 1e)	Managed and unmanaged grassland areas, including both intensive and extensive pastures.
Non-agricultural vegetation (Figure 1f)	Primary and secondary woody vegetation including paramos and closed-canopy forests
Water surface (Figure 1g)	Rivers, lakes, and artificial water surfaces (incl. reservoirs, ponds, and aquaculture).
Others (Figure 1h)	All other land use or land cover types, including built-up and open bare soil areas.
Other temporary crops (Figure 1i)	Areas with temporary crops including maize, soybean, and sugar cane in the preparation or harvest stage.
Other permanent crops (Figure 1j)	Coffee, cocoa, banana, and other permanent crop plantations include mixed systems with woody vegetation
Wetlands (Figure 1k)	Ecosystem that has geomorphologic conditions that cause saturation of the water in the surface.

To collect training points, we generated suitability maps for each class using secondary maps on forest cover (Hansen et al. 2013), water surface (Pickens et al. 2020), cropland (Potapov et al. 2021), palm oil (Descals et al. 2021), wetlands (SIAC 2020a), paramos (SIAC 2020b), and multiple LULC classes (González-González et al. 2022; IDEAM 2023) (Figs. S2 to S9 in the supplementary material). We further considered physical variables, such as the distance to water bodies, limiting suitable rice areas to regions of no more than 50 km distance, and slope in the case of grasslands used as pastures, which we assumed occur on slopes not exceeding 40%. (Zuluaga et al. 2021). We drew a stratified random sample (total n = 3,121) based on these maps and interpreted every point visually using very-high-resolution (VHR) images in Google Earth (GE). Uncertain or incorrect samples were removed, resulting in 1,623 labeled observations for 2020, 1,094 for 2015, and 404 for 2009. The suitability maps were used exclusively for model training, not for validation.

We extracted the STMs at our training locations and trained an initial Random Forest (RF) model in GEE using all 139 features (with 1,000 trees and trying $\sqrt{n_{\mathit{\text{features}}}}\mathrm{ }$ variables at each split) (Belgiu and Drăguţ 2016; Probst et al. 2019). The RF model was used to predict LULC classes across the study regions due to its documented performance in numerous remote sensing studies (Belgiu and Drăguţ, 2016; Rodriguez-Galiano et al., 2012). Random Forest classifiers are robust in high dimensional and multicollinear feature spaces, their parametrization is straightforward, and they allow estimation of uncertainty at the pixel-level, which enabled the use of active learning for improving model performance in our study (Belgiu and Drăguţ, 2016). We iteratively added points in selected locations (e.g. misclassified areas) manually and complemented the training sample with two runs of an active learning routine. Active learning helps to add highly informative training samples based on the calculation of model-based uncertainty (Stumpf et al. 2014). Following the routine described in Rufin et al. (2022a), we calculated probability margins derived from RF class probabilities as sampling strata to collect uncertain locations for training. Probability margins are calculated at the pixel level as the difference between the maximum class probability (i.e. class probability of the predicted class) and the second highest class probability. As such, large margins express high confidence scores, and low margins point to low model confidence. For each iteration, we used the class-wise 25th percentile of the distribution of probability margins in the study region to sample additional training points (20 points for each class) that were labeled using GE VHR imagery and PlanetScope mosaics provided through Norway’s International Climate and Forest Initiative (NICFI) data program (Planet Labs Inc 2020).

2.5. Area-adjusted accuracy assessment and area estimation

We followed the good practice recommendations for independent area-adjusted accuracy assessment based on a stratified random sample (Olofsson et al. 2014). We calculated the required sample size targeting a standard error of the overall accuracy (OA) of 1%, assuming user’s accuracies (UA) of 0.90 for the dominant classes, non-agricultural vegetation, and grasslands, and 0.85 for the other classes (Cochran 1977), yielding 937 samples. As we expected that generating reference labels for the classes of interest may be challenging given the limited availability of VHR imagery, we sampled a larger number of points (n = 2,559) for all years. We assigned at least 50 points per class and distributed the remaining samples according to class proportions, resulting in a higher number of samples for non-agricultural vegetation (n = 454), grasslands (n = 195), and others (n = 99).

All samples were interpreted using VHR imagery in GE and NICFI PlanetScope mosaics for the maps of 2015 and 2020. Samples were labeled independently by one experienced interpreter noting levels of confidence for the interpretation of each point for a second check by another experienced interpreter. Points flagged with high uncertainty were reexamined and labeled or deleted if no label could be assigned with certainty. We complemented the sample in four rounds of stratified random sampling to reach a sufficient number of validation samples. In total, we collected 1,396 for 2020, 1,392 for 2015, and 494 labeled samples for 2009. The relatively low number of labeled samples for the year 2009 was due to the scarcity of VHR imagery in parts of the study region, which increased standard errors in the estimation of accuracy and area. We generated a confusion matrix and derived area-adjusted OA, class-wise UA and PA, 95% confidence intervals, and error-adjusted area estimates from the reference data, using the mapac package for R v0.11 (Pflugmacher et al. 2019) using the estimators described in Stehman (2014).

2.6. Quantifying LULCC

We quantified the overall net LULCC based on the error-adjusted area estimates and conducted a spatial analysis in different subsets of the study region to investigate change trajectories. First, we considered all the changes that occur just in the area flooded by the reservoir of the dams to capture the immediate (direct) effects of reservoir flooding. Second, we considered all the municipalities that appear on the license as AOI (Figure 1, red area). Third, we analyzed LULCC trajectories in all municipalities in the upstream and downstream sub-basins of both dams, excluding the AOI (Figure 1, blue and brown area).

3. Results

3.1. Accuracy assessment and area estimates

The area-adjusted UA, PA, and OA with corresponding standard errors were calculated for each year and each class (Tables 3–5). The maps achieved satisfactory OA in 2020 (80.7%, ± 2.0%) and 2015 (80.0% ± 2.1%) and a medium OA (71.6% ± 3.9%) in 2009. UA ranging above 70% and partly exceeding 90% were registered for rice, palm oil, grasslands, non-agricultural vegetation, others, wetlands, and water surfaces. Producer’s accuracies higher than 70% and partly above 90% were registered for rice, palm oil, non-agricultural vegetation, water surface, and wetlands. On the contrary, temporary and permanent crops had low UA and PA across the years, indicating high omission errors for both classes.

Table 3. Confusion matrix for 2020 populated with probability scores, including UA and PA reported with standard errors.

Map class	Reference class (area proportion)									Accuracy
	2	3	4	5	6	7	8	9	10	Producer’s	User’s	F-score
Rice | 2	0.052	–	0.003	0.004	–	0.001	0.004	–	0.001	96.00 ± 2.22	79.12 ± 4.28	86.75 ± 2.73
Palm oil | 3	–	0.042	0.004	0.011	–	–	–	–	0.001	85.29 ± 4.26	72.50 ± 5.02	78.38 ± 3.44
Grassland | 4	–	0.001	0.099	0.007	–	0.004	0.003	–	0.003	56.79 ± 2.58	84.66 ± 2.77	67.98 ± 2.05
Non-agricultural vegetation | 5	0.001	0.006	0.043	0.409	0.001	0.009	0.012	0.019	0.002	91.21 ± 1.04	81.46 ± 1.46	86.06 ± 0.94
Water surface | 6	–	–	–	0.002	0.042	0.001	–	–	0.002	96.72 ± 2.25	88.06 ± 3.99	92.19 ± 2.41
Other | 7	–	–	0.005	0.006	–	0.064	–	–	–	78.26 ± 3.53	85.71 ± 3.43	81.82 ± 2.48
Other temporary crops | 8	0.001	–	0.012	–	–	0.001	0.029	–	0.001	59.42 ± 5.58	64.06 ± 6.03	61.65 ± 4.10
Other permanent crops | 9	–	–	0.005	0.009	–	–	0.001	0.021	–	51.79 ± 5.85	58.00 ± 7.04	54.72 ± 4.53
Wetlands | 10	–	–	0.003	0.001	0.001	0.001	–	–	0.049	81.93 ± 4.00	90.67 ± 3.38	86.08 ± 2.68
Overall accuracy = 80.66 ± 1.03

Table 4. Confusion matrix for 2015 populated with probability scores, including UA and PA reported with standard errors.

Map class	Reference class (area proportion)									Accuracy
	2	3	4	5	6	7	8	9	10	Producer’s	User’s	F-score
Rice | 2	0.038	–	0.001	0.001	–	–	0.002	–	0.001	96.36 ± 2.54	88.33 ± 4.15	92.17 ± 2.54
Palm oil | 3	–	0.024	–	0.003	–	–	–	–	–	70.83 ± 6.60	89.47 ± 5.01	79.07 ± 4.55
Grassland | 4	0.001	0.003	0.117	0.010	–	0.004	0.008	0.005	0.007	68.78 ± 2.82	75.12 ± 2.84	71.81 ± 2.01
Non-agricultural vegetation | 5	–	0.007	0.029	0.434	0.001	0.006	0.014	0.034	0.005	92.35 ± 1.01	81.73 ± 1.40	86.72 ± 0.91
Water surface | 6	–	–	–	0.002	0.032	0.001	–	–	0.005	91.67 ± 4.01	78.57 ± 5.52	84.62 ± 3.63
Other | 7	–	–	0.007	0.009	–	0.070	–	–	0.001	82.91 ± 3.28	80.83 ± 3.57	81.86 ± 2.43
Other temporary crops | 8	–	–	0.008	0.001	–	0.001	0.020	–	–	45.16 ± 6.08	66.67 ± 7.22	53.85 ± 4.92
Other permanent crops | 9	–	–	0.004	0.009	–	0.001	–	0.026	–	40.00 ± 4.88	63.16 ± 6.27	48.98 ± 4.12
Wetlands | 10	–	–	0.004	0.001	0.002	–	–	–	0.039	66.67 ± 5.10	85.71 ± 4.36	75.00 ± 3.63
Overall accuracy = 79.96 ± 1.06

Table 5. Confusion matrix for 2009 populated with probability scores, including UA and PA reported with standard errors.

Map class	Reference class (area proportion)									Accuracy
	2	3	4	5	6	7	8	9	10	Producer’s	User’s	F-score
Rice | 2	0.073	–	0.002	–	–	–	0.004	–	–	97.30 ± 2.70	92.31 ± 4.30	94.74 ± 2.60
Palm oil | 3	–	0.026	–	0.008	–	–	–	0.002	–	72.22 ± 10.67	72.22 ± 10.62	72.22 ± 7.53
Grassland | 4	–	0.006	0.136	0.022	–	0.008	0.004	0.006	0.004	60.91 ± 4.52	72.83 ± 4.52	66.34 ± 3.27
Non-agricultural vegetation | 5	–	0.002	0.040	0.273	–	0.008	0.022	0.047	0.004	81.82 ± 2.83	68.88 ± 3.28	74.79 ± 2.27
Water surface | 6	–	–	–	–	0.038	0.002	–	–	0.006	95.00 ± 4.81	82.61 ± 8.00	88.37 ± 5.03
Other | 7	0.002	–	0.010	0.024	0.002	0.091	0.004	–	–	83.33 ± 4.88	68.18 ± 5.74	75.00 ± 4.00
Other temporary crops | 8	–	–	0.016	–	–	–	0.018	–	–	34.62 ± 9.19	52.94 ± 12.12	41.86 ± 7.72
Other permanent crops | 9	–	–	0.004	–	–	–	–	0.010	–	15.62 ± 6.20	71.43 ± 17.19	25.64 ± 8.43
Wetlands | 10	–	0.002	0.014	0.006	–	–	–	–	0.051	78.12 ± 7.40	69.44 ± 7.44	73.53 ± 5.30
Overall accuracy = 71.66 ± 2.01

The confusion matrices reveal the dominant error types for each map (Tables 3–5). We observed high rates of confusion between permanent crops and non-agricultural vegetation, causing the above-mentioned omission errors for permanent crops. The dominant error types were similar across the three periods, except for the confusion between temporary crops and non-agricultural vegetation, which was particularly high in 2009, translating into the low producer’s accuracy for temporary crops.

Our error-adjusted area estimates reveal the extent of the target classes across the entire study region. Between the years the LULC classes changed dynamically but considering the whole period (2009–2020) the net changes were the decrease of wetlands, rice, grassland, temporary crops, and permanent crops, and an increase of water surface, non-agricultural vegetation, and palm oil (Figure 5).

Figure 5. Comparison of error-adjusted (filled bars) and map area estimates (black outlines). Error bars represent ± 1 standard error.

3.2. Spatial analysis of LULCC

The final LULC maps show the respective land cover according to our class catalog for the years 2009, 2015, and 2020 (Figure 6).

Figure 6. LULC maps for Hidrosogamoso (top) and El Quimbo (bottom) for 2009, 2015 and 2020.

According to our maps, the flooded area replaced by water surface in 2020 covered a total of 5,078.03 ha for Hidrosogamoso and 5,820.86 ha in the El Quimbo study region. For Hidrosogamoso, the flooded area dominantly replaced non-agricultural vegetation (77,9%) and grassland (14,3%) (Figure 7). The classes with the largest losses in the El Quimbo study region were rice (10,6%), grassland (25,9%), non-agricultural vegetation (47,7%), other (8,0%), and temporary crops (7,1%).

Figure 7. Alluvial plots showing the LULC trajectories for the reservoir areas hidrosogamoso (left) and El Quimbo (right) in 2020 for both study regions. The colors of the connecting flows correspond to the class in 2009.

LULCC in the AOI covered 0.7 Mha, or 14% of the Hidrosogamoso study area, and 0.2 Mha, or 6,3%, of the El Quimbo study region. Beyond the AOI, the changes in the Hidrosogamoso study area covered 4,2 Mha, and 2,7 Mha for the El Quimbo. The distribution of each LULC per year can be observed in Figure 8.

Figure 8. Alluvial plots showing the LULC trajectories of all changing pixels for (A) inside and (B) outside the AOI for the hidrosogamoso (left) and El quimbo (right) study region. Change trajectories below 10,000 pixels are omitted due to visualization purposes. The colors of the connecting flows correspond to the class in 2009.

In Hidrosogamoso, oil palm and grasslands suffered a loss inside and outside the AOI in both periods. Water surface decreased inside and outside the AOI in the first period, and increased in the second period. The other class decreased in both periods and spaces. Temporary and permanent crops increase and then decrease inside and outside the AI. Wetlands lost area in all regions and periods, despite in 2009–2015 in the AOI (Table 6). For El Quimbo, the rice areas suffered a decrease in the two periods and areas (inside and outside the AOI). The same happened for the other class. Grassland suffered decreases in all spaces and periods, despite 2009–2015. Temporary and permanent crops increase and then decrease for the influence area, and outside the influence area the first one decreases and then increases, and the second increases in both periods.

Table 6. Map-based LULC % change between the maps inside and outside the influence area.

	El-Quimbo				Hidrosogamoso
	Inside AOI		Outside AOI		Inside AOI		Outside AOI
Class	2009-2015	2015-2020	2009-2015	2015-2020	2009-2015	2015-2020	2009-2015	2015-2020
Rice	−0.8%	−0.2%	−0.1%	+0.0%	−0.1%	+0.1%	−0.1%	+0.0%
Oil Palm	+0.0%	+0.0%	+0.0%	+0.0%	−1.4%	+0.8%	−0.9%	+0.4%
Grassland	−4.7%	−4.2%	+0.8%	−3.1%	+1.9%	−2.5%	+1.7%	−2.0%
Non-agric. veg.	+1.5%	+9.1%	+0.2%	+4.3%	+5.8%	+2.2%	+2.1%	+3.0%
Water surface	+2.0%	+0.4%	−0.1%	0.1%	−0.5%	+0.1%	−0.3%	+0.1%
Other	−4.7%	−1.4%	−4.6%	−1.2%	−0.6%	−0.8%	−0.5%	−1.2%
Other temp crop	+0.3%	−1.5%	+0.5%	−0.6%	+0.1%	−0.1%	+0.1%	−0.1%
Other perm crop	+6.0%	−2.2%	+3.7%	0.5%	+0.2%	−0.1%	+0.5%	−0.2%
Wetlands	+0.1%	+0.1%	−0.3%	0.1%	−5.4%	+0.3%	−2.7%	−0.1%

3.3. Spatial analysis of cash crops

According to our maps, agricultural areas (i.e. the combination of temporary crops, permanent crops, rice, palm oil, and grassland) occupied 25% of the area in El Quimbo and 28% in the case of Hidrosogamoso region in 2020. This percentage increased between 2009–2015, mainly in the El Quimbo region, and maintained stable afterwards. For both cases, the predominant agricultural use was grassland (65% in El Quimbo and 81% in Hidrosogamoso), but note that this class includes unmanaged grassland as well. The grasslands were replaced over the three years mainly by palm oil in the case of the Hidrosogamoso and permanent crops in the case of El Quimbo. At the same time, our maps revealed a substantial decline in temporary crops. In the El Quimbo region, both downstream and upstream, this decrease reached 51,0%, between 2015 and 2020, and temporary crops were mainly replaced by water surfaces and grasslands (Figure 9A, B).

Figure 9. Examples of (A) areas flooded by the Hidrosogamoso reservoir, (B) areas flooded by the El Quimbo reservoir, (C) expansion of palm oil plantations at the expense of wetlands and grasslands, and D) palm oil areas affected by the disease. VHR imagery for 2009 in A (2009/12/05), B (2013/02/15), and C (2009/10/09) provided by ©maxar 2023 in Google Earth pro TM, VHR imagery for 2015 and 2020 in A (2016/06-11, 2020/03), B (2015/12, 2020/08), C (2016/12-2017/05, 2020/06-08), and D (2017/01/25, 2020/01) provided by Planet Labs Inc. ©2020.

The flooding of the Hidrosogamoso reservoir did not impact cash crops substantially, with most flooded agricultural lands being grasslands, and few temporary crop fields (Figure 9A). Palm oil covered the largest share of the agricultural area in the Hidrosogamoso region comprising 69.0% of the agricultural area in 2020. The palm oil area decreased between 2009 and 2015 but increased between 2015 and 2020. Palm oil has mainly replaced wetlands and grasslands (Figure 9C), with the former losing a substantial portion of their area between 2009 and 2015. On the other hand, palm oil has been replaced by non-agricultural vegetation (Figure 9D). Regarding cash crops in the El Quimbo region, rice and permanent crops dominated, occupying 37.2% and 50.3% respectively of the agricultural area in 2020. Rice areas decreased throughout the analysis period, particularly in the areas flooded by the reservoir (Figure 9B).

4. Discussion

4.1. Dam construction and land use/cover changes

We assessed LULC maps for three different periods to capture changes at different phases (construction, operation, present) and regions (flooded areas, licensed AOI, and beyond). We analyzed the changes in different spatial subsets of the region, to assess the changes beyond the areas of direct influence, and at the basin level.

The predominant classes in both study areas were non-agricultural vegetation and pastures. The presence of non-agricultural vegetation in much of both areas was also a pattern highlighted by González-González et al. (2022) in their mapping of the Andes region in 2018. Tendencies for forest regrowth or secondary vegetation succession were also found by Rubiano et al. (2017) in other parts of the Magdalena River basin. In this study, secondary vegetation mainly replaced marginal agricultural areas, while clearing dynamics were observed close to roads and urban areas. Pastures represented more than 70% of agricultural land uses, further reflecting the relevance of this extensive land use across Colombia (Etter et al. 2006; Zuluaga et al. 2021). In El Quimbo, reservoir flooding strongly affected agricultural land uses, contrary to Hidrosogamoso where non-agricultural vegetation occupied more than 75% of the flooded area; a pattern which is also confirmed by the license documents (ENEL 2008). The El Quimbo dam was constructed only 15 km upstream of the Betania dam, which was built as a multi-purpose dam in the 1960s. The joint use for irrigation and hydropower may have fostered the development of cash crops in this region as compared to the area surrounding Hidrosogamoso, where non-agricultural vegetation dominated, demonstrating how historic infrastructure expansion influences contemporary changing dynamics (De Andrade et al., 2021; De Oliveira Serrão et al., 2020).

Reservoir operation is linked to maintaining reservoir water levels according to the national energy demand. Thus, floodgates are closed or opened at different time periods and across seasons, affecting the river discharge and consequently water availability and accessibility downstream of the dams. As a consequence, we reason that the observed changes in wetlands and agricultural uses may relate to changes in the river flow controlled by the floodgates, providing additional insights into the interplay between dam commissioning, wetland conditions (Angarita et al. 2018) and downstream livelihoods (Baird et al. 2021). For Hidrosogamoso, wetlands areas decreased drastically (42,6%) between 2009 and 2015. This class was mostly replaced by grasslands downstream of the dams and in the influence area. Grassland expansion has been identified as an additional driver of wetland loss (Patino and Estupinan-Suarez 2016; Ricaurte et al. 2017), especially since this conversion is economically encouraged at the department level (Castiblanco et al. 2015; Ocampo-Peñuela et al. 2018).

For Hidrosogamoso, in 2015, a reservoir planning tool regulating competing land use interests around the reservoirs was published (‘Embalse Topocoro: Plan de Ordenamiento Central Hidroelétrica Sogamoso’). The document was elaborate by Fundación Humedales and ISAGEN, with the participation of social organizations. This instrument was required by Autoridad Nacional de Licencias Ambientales-ANLA (Resolution 1497/2009). LULC maps, as produced in this study, provide means to directly support management in the absence of planning instruments for El Quimbo, and enable modelling of related social, cultural, and economic drivers and consequences.

4.2. Cash crops, natural resources, and drivers of change

Colombia is the largest palm-oil producer in Latin America and the fourth-largest in the world (USDA 2023a). Palm oil dominates agricultural production in the Hidrosogamoso region, where oil palm cultivation was affected by a disease (pudrición del cogollo) starting in 2007 until 2013 (Rincón-Romero et al. 2021). Affected plantations were abandoned, until, in 2016, the disease could be controlled following technological developments and breeding of resistant species (Cabrera et al. 2020). The activation of formerly abandoned plantations has been captured by our maps with oil palm areas in 2009 changing to non-agricultural vegetation in 2015, followed by a partial change back to palm oil in 2020 (Figure 9D). Moreover, we identified municipalities with an increase in palm oil in the AOI, such as San Vicente de Chucuri, possibly related to the drastic price increase of palm oil in recent years (FEDEPALMA 2022). Currently, even small farmers have started to invest in palm oil. The creation and expansion of FUNDEPALMA, a federation dedicated only to smallholders, attest to these dynamics (Cabrera et al. 2020). These new oil palm plantations are established in proximity to rivers and are partly replacing wetlands, which is a new change trajectory as compared to the previous epoch, where oil palm mostly replaced grassland (Castiblanco et al. 2015).

Rice is among the most important cash crops in the El Quimbo region, occupying 34% of the agricultural area in 2020, which confirms the importance of Huila to national rice production (Agronet 2021). The areas concentrated near the main rivers and the dams remained stable during the three years. The retraction of rice areas evidenced in the maps occurred in the reservoir areas and areas at longer distances from water bodies, which are more vulnerable to changes in terms of water access. In 2010/2011, 23,759 ha of rice fields were impacted by La Nina, with a total loss of US$34.3 million, and in mid-2014, 33,938 ha were affected by El Niño-induced drought that led to crop losses valued at US$10.2 million (FEDEARROZ 2017). Recently, a retraction of rice areas has been documented in the Colombian Andes, which aligns with our findings (Castro-Llanos et al. 2019; Parra-Peña et al. 2022). Some of the remote rice areas were converted into water surfaces, which were visually identified as aquaculture areas. The Betania reservoir downstream of El Quimbo has been exploited by fish farming, following the initiatives by the departmental government to stimulate aquaculture production (Gobernación Huila 2006).

4.3. Uncertainties and limitations

Despite the data scarcity and dynamic nature of the Magdalena River basin, using multi-year STMs from Landsat images proved to be valuable for mapping LULC, even considering a complex class catalog. Our unbiased area estimates can be used for robust quantification of overall change dynamics in the region, while our maps are timely and due to the enhanced thematic detail in the agriculture domain, allow for more nuanced analyses of LULCC processes in the water-energy-food nexus, as compared to existing global and regional LULC products. Despite these advances, we identified remaining limitations, relating to limitations in data availability, the Landsat spatial resolution, and our multi-year aggregation strategy, which may have compromised the quality of our maps in certain areas. For specific periods and regions, we noted a low observation count and a seasonally clustered distribution of cloud-free observations, which likely reduced the distinguishability of the class-specific spectral-temporal behavior, mostly outside of the dams’ AOI.

Our area-adjusted accuracy assessment allowed for identifying the key error types, which in our case were commission errors of non-agricultural vegetation at the expense of temporary and permanent crops. The main permanent crop in the study region is coffee, which is mainly planted in oldgrown plantations (USDA 2023b), in agroforestry systems (Gomes et al. 2020), often in conjunction with other permanent crops such as banana or avocado (see Figure 1), which makes it appear similar to non-agricultural vegetation in the spectral-temporal feature space. Similarly, temporary crops have a high intra-class variability, while at the same time being spatially fragmented, and thus difficult to identify at Landsat spatial resolution. The limited data density required a multi-year aggregation strategy, which may have caused misclassifications for dynamic LULC types, including cropland-grassland rotations. On the contrary, the multi-year aggregation strategy is unlikely to affect rather stable LULC classes, such as rice. Such land uses require substantial land preparation and are often bound to the availability of surface water from rivers and streams for irrigation (Ochoa-Brito et al., 2023).

5. Conclusion

This study demonstrated that even in a data-scarce context like in Colombia, regional LULC maps derived from Landsat time series revealed changes across large regions and decadal time frames at high thematic detail. Integrating existing LULC products and iterative active learning proved to be a valuable tool to enhance training data collection, especially for complex classes in the agricultural domain. Our study provides further evidence that hydropower dams can directly affect LULC in different phases (construction and operation). Indirect effects can occur at spatial scales beyond the direct areas of influence, however, linking observed LULCC to the construction or operation of dams remains non-trivial, considering the interactions between multiple drivers of change (climate, social, cultural, economic). We thus argue that direct and indirect LULCC impacts need to be assessed and explicitly internalized in the dam’s environmental impact assessment and monitoring to be aligned with other territorial and policy instruments in order to decrease possible environmental and social inequalities.

Acknowledgements

This work was supported by the project ‘Water Security for Whom’, a collaboration between Humboldt-Universität zu Berlin (Germany), Pontificia Universidad Javeriana (Colombia), and Universidade Federal de Minas Gerais (Brazil). The project is funded by VolkswagenStiftung as part of the initiative ‘Global Issues – Integrating Different Perspectives on Social Inequality’ (project n. 96955). The authors would like to thank the ‘Water Security for Whom’ team, and the agricultural federations Federación Nacional de Arroceros (FEDEARROZ), Federación Nacional de Cultivadores de Palma de Aceite (FEDEPALMA), Federación Nacional de Cacaoteros (FEDECACAO) and Federación Nacional de Cafeteros de Colombia for meetings and checkings on the LULC maps during the fieldwork in Huila and Santander departments in Colombia in 2023. PR acknowledges funding by the F.R.S.-FNRS, grant no. T.0154.21.

Disclosure statement

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

Additional information

Funding

This work was supported by the project ‘Water Security for Whom’, a collaboration between Humboldt-Universität zu Berlin (Germany), Pontificia Universidad Javeriana (Colombia), and Universidade Federal de Minas Gerais (Brazil). The project is funded by VolkswagenStiftung as part of the initiative ‘Global Issues – Integrating Different Perspectives on Social Inequality’ (project n. 96955). The authors would like to thank the ‘Water Security for Whom’ team, and the agricultural federations Federación Nacional de Arroceros (FEDEARROZ), Federación Nacional de Cultivadores de Palma de Aceite (FEDEPALMA), Federación Nacional de Cacaoteros (FEDECACAO) and Federación Nacional de Cafeteros de Colombia for meetings and checkings on the LULC maps during the fieldwork in Huila and Santander departments in Colombia in 2023. PR acknowledges funding by the F.R.S.-FNRS, grant no. T.0154.21.

Notes on contributors

Caroline Salomão

Caroline Salomão contributed to Conceptualization, Methodology, Data curation, Software, Validation, Formal analysis, Investigation, Writing- Original draft preparation. Jonas Alsleben contributed to Methodology, Data curation, Software, Validation, Formal analysis, Writing-Review & Editing. Philippe Rufin contributed to Methodology, Writing-Review & Editing. Patrick Hostert contributed to Conceptualization, Writing-Review & Editing, Supervision.

Jonas Alsleben

Caroline Salomão contributed to Conceptualization, Methodology, Data curation, Software, Validation, Formal analysis, Investigation, Writing- Original draft preparation. Jonas Alsleben contributed to Methodology, Data curation, Software, Validation, Formal analysis, Writing-Review & Editing. Philippe Rufin contributed to Methodology, Writing-Review & Editing. Patrick Hostert contributed to Conceptualization, Writing-Review & Editing, Supervision.

Philippe Rufin

Caroline Salomão contributed to Conceptualization, Methodology, Data curation, Software, Validation, Formal analysis, Investigation, Writing- Original draft preparation. Jonas Alsleben contributed to Methodology, Data curation, Software, Validation, Formal analysis, Writing-Review & Editing. Philippe Rufin contributed to Methodology, Writing-Review & Editing. Patrick Hostert contributed to Conceptualization, Writing-Review & Editing, Supervision.

Patrick Hostert

Caroline Salomão contributed to Conceptualization, Methodology, Data curation, Software, Validation, Formal analysis, Investigation, Writing- Original draft preparation. Jonas Alsleben contributed to Methodology, Data curation, Software, Validation, Formal analysis, Writing-Review & Editing. Philippe Rufin contributed to Methodology, Writing-Review & Editing. Patrick Hostert contributed to Conceptualization, Writing-Review & Editing, Supervision.