Characterizing and Mapping the Grassland Vegetation of the Colombian Orinoquia

1. Introduction

Landscape-scale studies have shown that the distribution of plant communities is correlated with soil nutrients influenced by climate and geomorphology [1,2]. Bioclimatic indices, derived from quantitative temperature and precipitation data, are valuable tools for identifying thresholds associated with climatic events that induce hydric and thermal stress, thereby affecting vegetation patterns [3, 4]. Likewise, geomorphological factors are key determinants in vegetation modeling [5, 6, 7, 8]. Together with mesoclimatic patterns, these factors constitute the main predictors of the geographical distribution of vegetation along environmental gradients [9, 10, 11, 12, 13]. Mechanistic models that integrate spectral data with climatic variables such as temperature, solar radiation, and humidity have proven effective in the spatial characterization of vegetation [14]. This approach enables a better understanding of the ecogeographical ranges of plant communities, their responses to underlying environmental conditions, historical land-use changes, and intrinsic attributes such as floristic composition, characteristic species dominance, and autoecological traits, among others.

Grasslands occupy approximately 40% of the Earth's surface, making them one of the most extensive vegetation formations globally [15]. Their environmental characteristics and biotic composition provide key ecological, economic, and cultural services [16]. However, in recent decades, significant degradation processes linked to human activities have been documented, resulting in drastic alterations in species richness, floristic composition, and the structure of original plant communities [17].

The types of vegetation present respond to variables operating at three levels: (i) environmental factors or magnitude ranges—such as climate and topography—where phytocoenoses may establish; (ii) the assemblage of plant populations present, depending on their physiological responses; and (iii) the structural and functional characteristics of vegetation types [8]. Environmental complexes also influence the ecogeographic distribution of biodiversity in three main ways: (i) by acting as limiting or regulatory factors on the ecophysiology of organisms; (ii) through disturbances that modify environmental conditions; and (iii) as sources of assimilable resources or components [18, 19].

In this context, environmental and bioclimatic gradients are considered direct drivers, as they influence species physiology and, consequently, community metabolism. C₃ and C₄ grasses, for instance, exhibit significant differences in their distributional patterns along ecological, climatic, edaphic, and hydric gradients, due to substantial divergences in their physiology, biochemistry, anatomy, ultrastructure, and environmental requirements [20, 21, 22]. Topographic gradients, although considered indirect factors, are nonetheless strongly correlated with vegetation distribution, despite having a lesser direct influence on some biota [23, 24, 25, 26]. As a result, accurately representing the geographical distribution of phytocoenoses remains a persistent challenge, owing to their high spatial heterogeneity at multiple scales and their differential responses to climatic, geomorphological, and anthropogenic factors [27, 28].

This article presents a characterization of the current geographical distribution of grassland vegetation in the Colombian Orinoquia, based on the integration of data from: (i) georeferenced phytosociological alliances obtained from field surveys; (ii) the environmental complex, defined by the geomorphological gradient and limiting bioclimatic factors; and (iii) multisensor imagery, composed of Landsat-8 bands from the visible and infrared portions of the electromagnetic spectrum, along with dual-polarization Sentinel-1 bands. The aim is to improve the interpretation of these vegetation types in support of the development and implementation of conservation and sustainable management strategies.

2. Materials and Methods

2.1. Study Area

The Colombian-Venezuelan Llanos represent the second-largest expanse of Neotropical savannas, covering approximately 532,000 km² (254,000 km² in Colombia and 278,000 km² in Venezuela) [29, 30]. This region is notable for its floristic diversity and ecological complexity. In Colombia, it spans the departments of Arauca, Casanare, Meta, and Vichada, extending from the foothills of the Eastern Cordillera to the Orinoco River, and bounded by the Guaviare and Arauca rivers.

Its geographical boundaries are defined as follows: to the north, the Arauca River and the Venezuelan border; to the east, the Orinoco River; to the south, the Guaviare and Guayabero rivers; and to the west, the Guayabero River basin and the altitudinal limit of the Foothills (Piedemonte), which rises southward due to crustal flexure in the High Plain (Altillanura) [31].

From a geomorphological perspective, two main units are distinguished based on drainage and relief: the well-drained Orinoquia (east of the Meta River, in the departments of Meta and Vichada), characterized by dissected high plains, hills, and terraces; and the poorly drained Orinoquia (west of the Meta River, in Casanare and Arauca), marked by alluvial/eolic plains and alluvial fans [32]. Four physiographic units are recognized: the Foothills, Alluvial/Eolic Plains, High Plain, and La Macarena mountain range [32, 33, 34, 35].

2.2. Physiographic Setting.

The hydrographic zoning of Colombia [36] and the NASA SRTM v.3 digital elevation model (30 m resolution) were used within Google Earth Engine (GEE) to determine elevation ranges [37]. The Foothills zone [32] was delineated based on elevation thresholds by watershed: south of the Meta River (Guamal watershed), 575–350 m; north of the Meta, in the Meta and Pauto watersheds (675–200 m), Casanare (600–200 m), Cravo Sur (425–200 m), and Arauca (375–175 m). Its northern boundary was defined by the Arauca River and its southern boundary by the Ariari River watershed.

La Macarena was delimited by elevation thresholds of 575 m (Guayabero River) and 775 m (Guaviare River), and by the watersheds of the Ariari and Guayabero rivers. The Floodplain was identified to the west of the Foothills, bounded by the Arauca and Meta rivers and the eastern national border. The High Plain were defined as lying between the Foothills and La Macarena (to the west), the Meta River (north), the Orinoco River (east), and the Guaviare River (south) (Figure 1).

2.3. Grassland Phytosociological Framework [38].

The phytosociological classification comprises one class, three orders, and 18 alliances. The class Schizachyrio sanguinei–Trachypogonetea spicati is represented across all physiographic units of the region. It includes three orders: Axonopodo purpusi–Paspaletalia pectinati (comprising seven alliances), Schizachyrio brevifolii–Trachypogonetalia spicati (five alliances), and Axonopodo purpusi–Andropogonetalia virgati (six alliances). Their floristic composition, species richness, and ecological complexity span a wide range of environments—from piedmont foothills and dissected uplands in the high plains to aeolian flats and fluviolacustrine deposits within the alluvial plain.

2.4. Modelling Aproach. Data Sources. Mapping Framework

The adopted methodology follows an empirical–statistical modeling framework, based on the relationship between a multidimensional stack of predictor variables and grassland types as the response variable. The dataset consisted of 281 grassland relevés from the Colombian Orinoquia, compiled by Minorta-Cely et al., [38] and Rangel-Ch. et al., [39]. To model current distribution patterns, both unsupervised classification (K-means) and supervised classification (Random Forest) algorithms were applied. The former was used to define training areas, while the latter was used to predict spatial distribution. The modeling was supported by field data and integrated geomorphological variables, bioclimatic factors influencing vegetation distribution, and satellite imagery from Landsat-8 (optical) and Sentinel-1 (C-band SAR) sensors.

Bioclimatic variables were sourced from the WorldClim repository, which provides global interpolated temperature and precipitation data at 1 km resolution, derived from over 60,000 meteorological stations for the period 1970–2000 [40].

Figure 2 presents the methodological synopsis used to model the current distribution of vegetation. The input variables considered were: i) Geomorphology, which influences vegetation through factors such as solar exposure and the availability of water and nutrients; ii) Adjusted potential evapotranspiration, used as an indicator of thermal and hydric balance; iii) Annual ombrothermic index, which estimates water availability;

iv) Ombrothermic index for the driest quarter, which relates precipitation and temperature during that period; v) Thermicity index, which accounts for the influence of extreme temperatures; vi) C-band SAR data, which reflect physiognomic characteristics such as vegetation height and density; and vii) Optical sensor data, which capture the differential spectral response of grasslands in the visible and infrared ranges. The following section describes the cartographic layers used in the spatial modeling.

2.4. Data sources. Grassland Mapping

Geomorphology. Cartographic layer depicting the spatial distribution of regional-scale geomorphological units, derived from a model integrating key terrain descriptors: (i) HH/HV polarizations from an annual ALOS-PALSAR mosaic; (ii) NASA SRTM v.3 Digital Elevation Model (DEM); (iii) slope; (iv) aspect; (v) convexity; and (vi) roughness [32].

Adjusted Potential Evapotranspiration. Cartographic layer representing the maximum amount of water that can be evaporated from a vegetated soil under ideal conditions, adjusted by monthly insolation. Its calculation involved: (i) monthly mean temperature; (ii) annual heat index derived from monthly values, as a measure of accumulated thermal energy; (iii) the parameter a as a function of the heat index, recognized as an estimator of moisture retention in the air; and (iv) monthly solar radiation index calculated based on latitude and topographic characteristics derived from the NASA SRTM v.3 DEM, slope, and aspect.

Ombrothermic Index. Cartographic layer estimating annual water availability. The calculation incorporated: (i) annual positive precipitation and (ii) annual positive temperature.

Dry Quarter Ombrothermic Index. Cartographic layer estimating water availability during the driest quarter of the year, based on: (i) positive precipitation during the driest quarter, and (ii) positive temperature for the same period.

Thermicity Index. Cartographic layer quantifying the intensity of extreme temperatures. It was calculated using: (i) annual mean temperature; (ii) mean minimum temperature of the coldest month; and (iii) mean maximum temperature of the coldest month.

Sentinel-1 Mosaic. The Sentinel-1 SAR GRD collection provides dual-polarization C-band radar imagery, preprocessed, calibrated, and orthorectified using the Sentinel-1 Toolbox. This includes ground range detection, thermal noise removal, radiometric calibration, topographic correction using SRTM (30 m), and conversion to logarithmic decibel units [41]. Processing of 2022 imagery in Google Earth Engine (GEE) involved: (i) selection of the interferometric wide-swath mode as SAR source; (ii) filtering by ascending, descending, or both orbit directions depending on topography; (iii) filtering by VV and VH polarizations; (iv) selection of 10 m spatial resolution data; (v) clipping to the area of interest uploaded as an Asset in GEE; (vi) alternating selection of VV and VH bands; (vii) full-year temporal filtering; (viii) mosaic composition from both bands; and (ix) application of a 50 m circular focal mean filter to reduce speckle noise.

Landsat-8 Mosaic. This mosaic consists of 2022 images from the USGS Landsat 8 Surface Reflectance Tier 1 collection, which includes bands from the visible spectrum (RGB), near-infrared (NIR), and shortwave infrared (SWIR), useful for vegetation differentiation [42, 43, 44]. The collection [45] contains data atmospherically and thermally corrected using the OLI/TIRS sensor and LaSRC software. Cloud and shadow masking was applied using the QA band, and the mosaic was generated using the pixel-wise median. NDVI was calculated from red and NIR bands and incorporated into the model to reduce the effects of illumination, topographic heterogeneity, and soil reflectance [46]. The GEE processing steps included: (i) application of a cloud and shadow masking function using the QA band; (ii) calculation of the median of unmasked pixels; and (iii) NDVI computation and stacking with the remaining bands.

Image Processing and Classification in GEE. Google Earth Engine (GEE) allows for analysis without the need for explicit scale or projection definitions, as these are either automatically inferred or determined by the output data [47]. However, to avoid inconsistencies in multi-source mosaics, both projection and scale were defined manually. The procedure included: (i) reprojection to 30 m resolution based on native Landsat-8 data; and (ii) stacking of 12 bands into a multidimensional image: one geomorphological unit layer, four bioclimatic variables, two radar bands (VV + VH), and five optical bands (RGB, NIR, and NDVI).

Training Area Definition. Training areas were delineated using an unsupervised classification approach with the K-means algorithm, based on stacked image data and the spatial location of vegetation alliances identified during field surveys. K-means partitions the dataset into k clusters by assigning each observation to its nearest centroid [48]. The procedure included: (i) pixel sampling according to geometry, scale, and desired quantity; (ii) algorithm training with the number of vegetation alliances plus four general land cover types (forest, agricultural use, urban areas, and water bodies); (iii) classification of the multidimensional image into a single-band raster containing cluster identifiers; and (iv) export at a spatial resolution of 30 m.

The classified image was generalized in GIS software by reclassifying non-vegetated covers as “NoData” and aggregating pixels into minimum mapping units (MMUs) of 5 hectares. This involved: (i) thematic aggregation using a majority focal filter (2x2 window); (ii) clumping of contiguous pixels by class; and (iii) iterative elimination of areas smaller than 5 ha by reassigning them to the dominant neighboring class. The final raster was vectorized and spatially linked to vegetation survey points.

Random Forest Classifier Configuration. The configuration of the Random Forest classifier was based on the optimal number of decision trees and training variables. While increasing the number of trees generally improves model accuracy, it also raises computational costs [49]. The tuning process involved: (i) sample extraction from the multidimensional image and training polygons; (ii) division of the dataset into 70% for training and 30% for validation [50]; (iii) iterative testing with varying numbers of trees to optimize accuracy while avoiding overfitting [51, 52]; (iv) classifier configuration using the optimal tree count; (v) training with 70% of the data and classification using numerical attributes and stacked bands; (vi) validation with the remaining 30%, including confusion matrix computation and accuracy plots as a function of tree number. The final classification employed the optimal number of trees and the same set of training variables [53]. The result was exported as a raster at 30 m resolution, generalized to 5 ha MMUs, and vectorized. Accuracy metrics were computed, including the confusion matrix, overall accuracy, producer’s and user’s accuracies, and the Kappa coefficient.

2.5. Grassland Ecological Atributes

For each phytosociological alliance, two intrinsic components were considered: (i) the classification of photosynthetic pathways (C₃ or C₄) of its characteristic grasses, which was addressed from a foliar anatomical perspective, based on observations using optical microscopy (OM) and scanning electron microscopy (SEM). Analyses were performed on transverse sections of mature leaves and paradermal views of leaf fragments [21, 54, 55]. (ii) Dominance Ratio (DR_i), defined as the cardinal proportion of dominance based on the following premise: the ratio between the number of dominant taxa (i) and the total number of taxa comprising the alliance. Let the following sets be defined:

S: the total set of taxa present in a given alliance.

D ⊆ S: the subset of dominant taxa within the alliance.

Then: |S|: the total number of taxa i (i.e., the cardinality of set S).

|D|: the number of dominant taxa i (i.e., the cardinality of subset D).

Thus: DR_i = (|D| / |S|) * 100. Where DR_i represents the percentage of dominant taxa relative to the total taxa constituting the alliance, expressed as a percentage within the range [0, 100].

3. Results

The results indicate that 73.74% of the territory retains some form of natural vegetation. Of this proportion, 32.19% corresponds to forests, while 41.55% comprises 18 phytosociological alliances of natural grasslands. The most extensive alliances include Paspalo pectinati–Axonopodion aurei (6.02%), Axonopodo aurei–Trachypogonion spicati (4.37%), and Rhynchosporo barbatae–Axonopodion ancepitis (4.44%). In contrast, the alliances Sipaneo pratensis–Axonopodion purpusi, Rhynchosporo corymbosae–Schyzachyrion brevifolii, and Schizachyrio brevifolii–Tibouchinion asperae collectively account for less than 2% of the total area (Table 1).

3.1. Foothills Graslands

13.94% of the unit preserves natural vegetation, with 12.51% corresponding to forests and 1.43% to natural grassland formations (Figure 3). Five phytosociological alliances were recorded, notably Eleochario interstinctae–Rhynchosporion barbatae (1.05% cover) and Trichanthecio cyanescentis–Andropogonion virgati (0.33%), both located on accumulation terraces along the eastern margin of the physiographic unit, bordering the Orinoquia Floodplain.

The transformation of 86.06% of the natural vegetation in this unit is primarily due to human activities, particularly extensive cattle ranching and agro-industrial crops, especially African oil palm (Elaeis guineensis Jacq.) and rice [56]. The physiographic conditions—gentle slopes and sediment inputs from the Andes—have facilitated the expansion of agricultural and livestock activities over denudational and fluvial environments. In particular, agriculture has expanded through the establishment of oil palm plantations on sub-recent terraces, while mining, concentrated in Floodplain areas, is focused on fluvial gravel extraction. Areas with the highest anthropogenic impact, such as urban infrastructure and road networks, are concentrated on alluvial fans, taking advantage of their accessibility and proximity to agricultural lands [32].

3.2. La Macarena Grasslands

60.62% of the unit retains natural vegetation, with 55.22% corresponding to forest and the remainder to open vegetation. The alliance Paspalion carino-pectinati covers 3.88% of the area and is primarily distributed across dissected hill systems in the lower Ariari River basin, near the boundary with the High Plain (Figure 4). The remaining 39.38% of the area has been transformed, mainly due to extensive cattle ranching on low-fertility soils, which has driven deforestation and land concentration in areas with poor agricultural suitability. Other anthropogenic pressures contributing to this transformation include subsistence farming, illicit crops, agricultural frontier expansion, burning practices, and human settlements in high-risk zones—trends indicative of deficient territorial planning [57, 58]. These pressures particularly affect non-flooded forests and vulnerable hydrographic networks, which are preferred for agricultural use due to their well-drained soils and lower acidity [59].

3.3. Floodplain grasslands

66.14% of the unit maintains natural vegetation, of which 10.95% corresponds to forest and 55.20% to natural grasslands (Figure 5). Seven syntaxa were recorded, with Rhynchosporo barbatae–Axonopodion ancepitis being the most widespread (19.79%), primarily associated with aeolian geomorphological units located along the eastern margin. Other notable alliances include Eleochario interstinctae–Rhynchosporion barbatae, common on terraces and denudational plains in the central and southern sectors; Trichanthecio cyanescentis–Andropogonion virgati, found on terraces adjacent to eolic landscapes; and Steinchismo laxae–Andropogonion bicornis, located in complex fluvial environments in the northwestern region, where denudational plains converge with floodplains. Less representative grassland types include Rhynchosporo barbatae–Andropogonion virgati, Andropogono virgati–Axonopodion ancepitis, and Caperonio palustris–Leersion hexandrae, located on terraces and aeolian formations along the western margin of this unit.

Human activities have altered 33.86% of the natural vegetation in this unit. Agricultural zones, including subsistence crops and introduced pastures, are concentrated in low-lying areas (bajos/bajíos), floodplains (vegas), and seasonally flooded savannas (vegones) associated with the Lipa, Ele, and Cravo Norte rivers and their tributaries, where soils exhibit lower acidity and reduced aluminum content compared to those of the terraces [60].

3.4. High Plain Grasslands

86.45% of this unit retains natural vegetation, with 39.30% corresponding to forest and 47.15% to open vegetation (Figure 6). Eleven grassland alliances were identified, the most representative being Paspalo pectinati–Axonopodion aurei (9.94%), which is distributed across peneplains and gently undulating plains between the Meta, Muco, Tomo, and Bita rivers, as well as north of the Tuparro River. It is followed by Axonopodo aurei–Trachypogonion spicati, which occurs on dissected hill systems south of the Cumaral stream, along the eastern margin of the Manacacías River, and in the Upper Vichada region on collinated plains associated with the Tomo, Bita, and Tuparro rivers.

Hyptio confertae–Schizachyrion brevifolii has a broad but fragmented distribution, except in the southern portion of the Vichada River basin, where forest cover predominates. Paspalion carino–pectinati is common in the western sector, particularly in highly dissected hills between tributaries of the Manacacías and Guaviare rivers, and in the northern part of the upper Vichada basin. Xyrio savanensis–Coletaenion caricoides occurs in deeply dissected hills near the Yucao River, in the middle Vichada basin, and in denudational zones influenced by the Tomo and Bita rivers.

Rhynchosporo corymbosae–Schyzachyrion brevifolii is prevalent in hill systems near the Manacacías and Vichada rivers. Other alliances with lower representativity and scattered distributions include Hyptio brachiatae–Trachypogonion spicati, Sacciolepio angustissimatis–Rhynchanterion bracteatae, Schizachyrio brevifolii–Tibouchinion asperae, Sipaneo pratensis–Axonopodion purpussi, and Hyptio lantanifoliae–Ichthyotherion terminalis.

The remaining 13.55% of the area has been transformed by human activities, particularly in fluvial environments where chemical reduction processes decrease soil acidity and improve fertility, thus favoring the expansion of agricultural and livestock activities [61].

3.5. Ecogeographical Considerations

3.6. Dominance Ratio (DR_i)

The phytocoenoses are composed, on average, of 51 diagnostic species and eight characteristic species, meaning that approximately 20% of their species exhibit significant dominance values (DR_i). The ratio of Poaceae dominance relative to other characteristic plant families (DR_P) ranges from 16.7% (Hyptio lantanifoliae–Ichthyotherion terminalis) to 75% (Paspalo pectinati–Axonopodion aurei). Only Eleochario interstinctae–Rhynchosporion barbatae lacks characteristic Poaceae species. Table 2.

3.8. Photosynthetic Dominance

There is a clear dominance of Poaceae species with C₄ photosynthetic pathways. Rhynchosporo barbatae–Andropogonion virgati and Sipaneo pratensis–Axonopodion purpusi are particularly noteworthy, as they not only show a marked dominance of Poaceae among their characteristic species (50% and 60%, respectively), but all of these are C₄ species.

4. Discussion

Multisensor applications facilitate the development of novel technological approaches and overcome inherent technical limitations [62]. The integration of field data with satellite imagery helps bridge the gap between remote sensing and ecological research, promoting efficient methods for large-scale, continuous, and detailed biodiversity assessments [63, 64, 65, 66]. The integration of multisource and multitemporal data was crucial for discriminating vegetation types, as it provided complementary information on the geometry and texture of vegetation units, thereby reducing the uncertainty associated with structural and spectral similarities [67, 68, 69, 70].

4.1. On the Catographic Model

The spatial discrimination of grassland vegetation was enhanced through the integration of data from diverse sensors and thematic layers, which provided complementary information on the geometry and texture of alliances. This integration reduced uncertainty caused by structural and spectral similarities [67, 68, 69, 70]. The accuracy of classification depends on the quality of the training samples, which must adequately represent the spectral and landscape variability of the grasslands [71, 72]. Balanced samples across thematic classes increase classification accuracy by reducing commission and omission errors [46, 73]. The expansion of spectral and thematic resolution through multidimensional imagery improved overall accuracy and reduced classification variance [74, 75]. The use of the annual per-pixel median mosaic mitigated phenological effects related to seasonal alternation, thereby enhancing classification precision [72].

Random Forest models were used to classify stacked images into discrete units according to the pre-delimited grassland alliances. Optimal tree parameterization ranged from 90 (Foothills) to 160 (High Plain), yielding overall accuracies between 0.870 and 0.790. These values indicate a 79–87% match between classified pixels and reference categories. Confusion matrices were employed to evaluate classification accuracy by examining the agreement along the main diagonals, error distribution, and the presence of non-random pixel assignments. Maximum omission and commission errors were 25% and 31%, respectively, and were primarily associated with the High Plains region. The Kappa coefficient, which measures the agreement between automated and random classifications [76, 77], reached minimum values of 0.91, indicating a high probability of correct classification. This digital approach enhances consistency and reproducibility, avoiding the subjectivity inherent in visual classifications [64, 69].

Although K-means proved useful in defining representative training areas, future research could adopt stratified sampling schemes that select balanced training data based on updated thematic maps. This approach would improve the accurate inclusion of training areas without relying on their spatial extent [72, 73, 78]. Additionally, Random Forest could be employed to rank variables based on their discriminative contribution, facilitating the refinement or incorporation of new input variables [51, 72, 79].

4.2. Dominance Ratio (DR_i) and Photosynthetic Pathways of Poaceae

C₄ grasses are more diverse in tropical lowlands and tend to dominate open habitats, reaffirming their prominence in terms of both species richness and cover within phytocoenoses. This marked presence of C₄ species is directly related to their functional and physiological traits. In general, these species are characterized by a low atmospheric CO₂ compensation point, rapid growth rates, low rates of photorespiration and water loss per unit of biomass produced, a broad optimal temperature range, a high light saturation point, and a foliar Kranz-type anatomy. In contrast, C₃ grasses exhibit a high CO₂ compensation point, elevated photorespiration rates, an optimal temperature range between 10 and 25 °C, a substantially lower light saturation point than C₄ species, and non-Kranz-type leaf anatomy [54].

These traits help explain, at least in part, the dominance of C₄ species in various environments of the Orinoquia region and the limited presence of C₃ grasses. The low richness and ecological representation of C₃ grasses in the Orinoquia is notable, as the few species exhibiting this photosynthetic pathway in the region are almost exclusively found in shaded habitats (e.g., understories and ecotones between forests and open habitats, or aquatic environments). This distribution correlates with one of the characteristic limitations of C₃ species—their lower light saturation point.

As such, C₄ plants hold a competitive advantage over C₃ species in environments with high light intensity, elevated temperatures, high atmospheric O₂ concentrations, and low CO₂ concentrations. Conversely, the higher energy cost associated with the C₄ metabolic pathway, compared to the C₃ pathway, excludes C₄ species from shaded environments and reduces their competitive advantage in colder habitats. Therefore, the distribution of C₄ grasses is largely climate-driven and corresponds primarily to open habitats [20].

The absence of C₄ grasses in closed vegetation may be due to insufficient evolutionary time for this photosynthetic pathway to adapt to such environments. The C₄ pathway evolved relatively recently under warm, and generally arid, conditions. It has been found that C₄ photosynthesis developed in tropical-origin grasses during the Miocene, a period in which C₄ grasses also expanded through multiple evolutionary originssee references in [20, 54]. Additionally, van Els et al., [80] suggested that temporal estimates for the emergence of tropical grasslands dominated by C₄ grasses range from the early to the late Miocene, with the latter interval being particularly relevant for South America. Neotropical grasslands originated during the Paleogene; however, those dominated by C₃ species began to diversify by the mid-Cenozoic, whereas the expansion of C₄-dominated grasslands occurred later, during the late Neogene which coincides with the geomorphological development of the Orinoquia region.

5. Conclusions

In Colombia, biodiversity conservation is often guided by surrogate indicators, such as land cover maps. Although useful, these proxies can be problematic when their correspondence with the actual distribution of biodiversity is low [81, 82]. Detailed thematic mapping has proven more effective in addressing environmental conflicts and enhancing biodiversity management strategies [1, 64]. Vegetation maps serve as essential inputs for wildlife habitat studies [83], species and community distribution modeling [81, 84, 85], identification of priority areas for conservation [1], and the evaluation of threats and disturbances to biodiversity [86, 87]. They also support the sustainable management of forest resources [88] and the biophysical estimation of environmental supply. Furthermore, they inform ecological restoration and compensation processes by providing accurate information on dominant vegetation composition [28].

The current distribution of grassland vegetation in the Colombian Orinoquia was characterized using machine learning techniques applied to multidimensional imagery, through an empirical-statistical modeling approach based on the relationship between stacked image data (predictor variables) and grassland alliances (response variable). A total of 73.74% of the region still retains some form of natural vegetation. Of this, 32.19% corresponds to forest, while 41.55% comprises 18 natural grassland syntaxa. The main anthropogenic activities impacting natural grasslands include agricultural expansion, subsistence farming, illicit crop cultivation, and the introduction of forage pastures based on introduced or naturalized species.