Predictive Assessment of Climate Change Impact on Water Yield in the Meta River Basin, Colombia: An InVEST Model Application

: This paper presents a hydrological assessment of the 113,981 km 2 Meta River basin in Colombia using 13 global climate models to predict water yield for 2050 under two CMIP6 scenarios, SSP 4.5 and SSP 8.5. Despite mixed performance across subbasins, the model was notably effective in the upper Meta River subbasin. This study predicts an overall increase in the basin’s annual water yield due to increased precipitation, especially in ﬂatter regions. Under the SSP 4.5, the Meta River basin’s water ﬂow is expected to rise from 5141.6 m 3 /s to 6397.5 m 3 /s, and to 6101.5 m 3 /s under the SSP 8.5 scenario, marking 24% and 19% increases in water yield, respectively. Conversely, the upper Meta River subbasin may experience a slight decrease in water yield, while the upper Casanare River subbasin is predicted to see signiﬁcant increases. The South Cravo River subbasin, however, is expected to face a considerable decline in water yield, indicating potential water scarcity. This study represents a pioneering large-scale application of the InVEST–AWY model in Colombia using CMIP6 global climate models with an integrated approach to produce predictions of future water yields.


Introduction
Climate change is a global issue that has widespread impacts on various aspects of the environment, including water resources.The effects of climate change on water resources are multifaceted and can have significant consequences for freshwater ecosystems, biodiversity, and human societies [1-3].One of the key impacts of climate change on water resources is the alteration of precipitation patterns [2].Dai [4] highlights that climate change results in changes in precipitation patterns affecting water yield and leads to water scarcity or excess in different areas [1].Water yield refers to the total amount of water that is generated within a specific drainage basin or watershed over a given period of time according to Budyko [5].Freshwater ecosystems are particularly vulnerable to the effects of climate change on water yields.
Hydrological models are widely used to study the impact of climate change on the hydrologic cycle and the management of water resources systems [6].These models integrate various hydrological and ecological processes that allow for the simulation of water yield under different environmental and climatic conditions [7], providing valuable insights into the potential effects of climate change scenarios [8].However, many of these models require amounts of data often not available in developing countries, which are likely to be the most affected by climate change.
Therefore, there is a growing need for assessing water yield response in areas with low weather station densities, such as the tropics, through hydrological models with low data requirements.Despite limitations assessing model performance, several studies have shown that relatively simple models with limited data requirements can credibly simulate flows and water balance components in large, data-poor basins.Valencia et al. [9] used the annual water yield model from the INVEST package to quantify current water yields in specific subbasins in the Meta River basin, Colombia.The study [10] also demonstrated the application of a distributed hydrological model like MIKE-SHE in data-scarce watersheds to assess climate change impacts on water resources, using limited meteorological and hydrological data for model input.
These simple models have been used to analyze the effects of climate change on water yields in different regions.The IHACRES model, used by Kim et al. [11], was effectively applied in various catchments, accounting for nonstationary and seasonal effects in climate impact studies.Lastly, models requiring low data, as used by Bejagam et al. [12], demonstrated effective results from the InVEST annual water yield model in assessing the impacts of climate change on water yield in the Tungabhadra River basin, India.
A limited number of studies have analyzed the impact of climate change on Colombia's water resources at the basin scale.Various climate change scenarios predict the country's rainfall to increase, but ongoing complex land-use dynamics and data shortages create uncertainty about the impacts of climate change in certain basins.The Meta River basin has seen a significant expansion of agricultural activity and has been the center of many development plans in Colombia; however, hydrological and climatic data in the region are limited.
The primary objective of this research was to perform a future prediction of the annual water yield in the Meta River basin (Colombia) by 2050.Building on work [9] that revealed historical and current patterns of water yield, this assessment was conducted using the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model for annual water yield.This research emphasized analyzing water yield predictions under two specific Shared Socioeconomic Pathways (SSPs) of the Coupled Model Intercomparison Project Phase 6 (CMIP6)-4.5 and 8.5 scenarios.A feature of this study was the novel use of an ensemble approach combining 13 global climate models under these two CMIP6 scenarios.This methodology increased the robustness of the predictions, allowing for a wide assessment of future water availability in the Meta River basin.Based on foundational research previously carried out in this area by [9], this allowed for a comprehensive predictive assessment of large-scale water availability in the region, which is crucial for Colombia's future agricultural productivity.

Study Area
The boundary of the Meta River basin was defined by [9], with the aid of ArcSWAT tool and the Digital Elevation Model (DEM) having a 30 m resolution, obtained from the Global Multi-Resolution Topography (GMRT) dataset [13].This process outlined a region spanning 113,981 km 2 (Figure 1).In addition, the InVEST-AWY model's effectiveness was assessed in four specific subbasins of the river, chosen based on the presence of gauging stations.The purpose of this assessment was to identify the regions where the model produced the most accurate results.The Meta River, a left tributary of the Orinoco River, traverses several Colombian departments, including Meta, Casanare, Cundinamarca, Boyacá, Arauca, Vichada, and Bogotá.This region includes 31.7% of the country's flood-prone zones [14] and has a tropical climate and varied topography, including mostly flat and gently sloping lands, as well as mountainous regions in the west, including the eastern Cordillera and La Macarena range.Air temperature varies between 4 • C and 28 • C, with annual rainfall ranging from 1000 to 7000 mm [14,15].Climatic conditions within the region vary: the lowland areas have a monomodal climate with a single rainy season, contrasting with the bimodal climate of the mountains [16].
and gently sloping lands, as well as mountainous regions in the west, including the eastern Cordillera and La Macarena range.Air temperature varies between 4 °C and 28 °C, with annual rainfall ranging from 1000 to 7000 mm [14,15].Climatic conditions within the region vary: the lowland areas have a monomodal climate with a single rainy season, contrasting with the bimodal climate of the mountains [16].

Method
The methodology flowchart for this study, adapted from [9,12], is shown in Figure 2. The initial step involved preparing key spatial datasets such as precipitation, potential evapotranspiration, plant available water content, root restricting layer, land-use/landcover (LULC) with a biophysical table, and basin/subbasin delimitation maps.The sensitivity analysis carried out by [9] used a base Z (the Zhang coefficient) value of 30, which is the annual number of rainfall days divided by 5.In the study area, the average number of rainfall days is 177.The coefficient Z was adjusted in the range from 1 to 30 (Figure 3).In relation to the study region, our results showed that at Z = 15, water yield increased by 9%, while at Z = 1, it increased by 101%.

Method
The methodology flowchart for this study, adapted from [9,12], is shown in Figure 2. The initial step involved preparing key spatial datasets such as precipitation, potential evapotranspiration, plant available water content, root restricting layer, land-use/landcover (LULC) with a biophysical table, and basin/subbasin delimitation maps.
and gently sloping lands, as well as mountainous regions in the west, including the eastern Cordillera and La Macarena range.Air temperature varies between 4 °C and 28 °C, with annual rainfall ranging from 1000 to 7000 mm [14,15].Climatic conditions within the region vary: the lowland areas have a monomodal climate with a single rainy season, contrasting with the bimodal climate of the mountains [16].

Method
The methodology flowchart for this study, adapted from [9,12], is shown in Figure 2. The initial step involved preparing key spatial datasets such as precipitation, potential evapotranspiration, plant available water content, root restricting layer, land-use/landcover (LULC) with a biophysical table, and basin/subbasin delimitation maps.The sensitivity analysis carried out by [9] used a base Z (the Zhang coefficient) value of 30, which is the annual number of rainfall days divided by 5.In the study area, the average number of rainfall days is 177.The coefficient Z was adjusted in the range from 1 to 30 (Figure 3).In relation to the study region, our results showed that at Z = 15, water yield increased by 9%, while at Z = 1, it increased by 101%.The sensitivity analysis carried out by [9] used a base Z (the Zhang coefficient) value of 30, which is the annual number of rainfall days divided by 5.In the study area, the average number of rainfall days is 177.The coefficient Z was adjusted in the range from 1 to 30 (Figure 3).In relation to the study region, our results showed that at Z = 15, water yield increased by 9%, while at Z = 1, it increased by 101%.
In contrast, the sensitivity analysis for the crop evaporation factor, Kc, revealed that a 30% reduction in Kc resulted in a 10.7% decrease in water yield.Conversely, increasing Kc by 30% only resulted in an increase in water yield of only 2%.These results suggest that, in the InVEST-AWY model, the Kc value does not have a significant effect on water yield, unlike the Z value.In contrast, the sensitivity analysis for the crop evaporation factor, Kc, revealed that a 30% reduction in Kc resulted in a 10.7% decrease in water yield.Conversely, increasing Kc by 30% only resulted in an increase in water yield of only 2%.These results suggest that, in the InVEST-AWY model, the Kc value does not have a significant effect on water yield, unlike the Z value.

The InVEST-AWY Model
The InVEST Annual Water Yield model is designed to estimate how different parts of the landscape contribute to water availability.It offers valuable perspectives on the impact of changes in land-use and land-cover on annual surface water yield and hydropower generation, as highlighted in [17].The water yield module in the InVEST model is based on the annual average precipitation and the Budyko curve [5].The annual water yield for each pixel is determined by the Equation (1): where AET(x) is the annual average evapotranspiration, and P is annual precipitation for each pixel (x).For lands with known vegetation or land-use/cover types (LULC), the evapotranspiration fraction of the water balance is ( ) ( ) ; it is calculated as follows: where PET(x) represents the annual potential evapotranspiration; w(x) is a parameter characterizing soil-climate properties; the formulas for PET(x) and w(x) are as follows: where Kc(x) is the crop coefficient per pixel;  is the potential evapotranspiration per pixel; w(x) is a non-physical parameter that characterizes the natural climatic properties of the soil (Equation ( 4)); () is the water available to the plant; Z is the Zhang coefficient, which depends on the amount of precipitation per year.For other land-use and land-cover types like open water surfaces, urban areas, and wetlands, actual evaporation is directly calculated from the reference evaporation,  (), with an upper limit set by precipitation as described in Equation (5):

The InVEST-AWY Model
The InVEST Annual Water Yield model is designed to estimate how different parts of the landscape contribute to water availability.It offers valuable perspectives on the impact of changes in land-use and land-cover on annual surface water yield and hydropower generation, as highlighted in [17].The water yield module in the InVEST model is based on the annual average precipitation and the Budyko curve [5].The annual water yield for each pixel is determined by the Equation (1): where AET(x) is the annual average evapotranspiration, and P is annual precipitation for each pixel (x).For lands with known vegetation or land-use/cover types (LULC), the evapotranspiration fraction of the water balance is AET(x) P(x) ; it is calculated as follows: where PET(x) represents the annual potential evapotranspiration; w(x) is a parameter characterizing soil-climate properties; the formulas for PET(x) and w(x) are as follows: where Kc(x) is the crop coefficient per pixel; ET 0 is the potential evapotranspiration per pixel; w(x) is a non-physical parameter that characterizes the natural climatic properties of the soil (Equation ( 4)); AWC(x) is the water available to the plant; Z is the Zhang coefficient, which depends on the amount of precipitation per year.For other land-use and land-cover types like open water surfaces, urban areas, and wetlands, actual evaporation is directly calculated from the reference evaporation, ET 0 (x), with an upper limit set by precipitation as described in Equation (5): where ET 0 (x) is the reference evapotranspiration; Kc(lx) denotes the evaporation factor specific to each land-use and land-cover type.

Data Requirements
To calculate the annual water yield of the river with the InVEST-AWY model, it is essential to input various raster and shape format data (Table 1).These include the annual average rainfall, annual average reference evapotranspiration, land-use and land-cover data accompanied by a biophysical table, the depth of the layer that restricts root growth, the amount of water available for plant use, and detailed maps of the studied river basin and its subbasins.

Meteorological Data
Meteorological data were sourced from the IDEAM (Instituto de Hidrología, Meteorología y Estudios Ambientales) website [20], including figures on annual rainfall (Figure 4A), annual average water flow, and annual average maximum and minimum air temperatures.This study identified a total of 246 hydrometeorological stations in the upper Meta River subbasin that record air temperatures, precipitation, and water discharge.In contrast, the Yucao River subbasin had only one station, and the South Cravo River subbasin had four stations monitoring precipitation and air temperatures.In addition, there were 20 hydrometeorological stations in the upper Casanare River subbasin.Records from these in situ gauging stations date back to 1983.The annual potential evapotranspiration (Figure 4E) was estimated using the Hargreaves formula, which uses air temperature data from on-site stations and extraterrestrial solar radiation data derived from [21], using the R package environment.These data were spatially interpolated to a resolution of 1 km × 1 km.The potential evapotranspiration-PET (Eto)-was calculated as follows: where Tmax and Tmin are annual average maximum and minimum air temperatures ( • C), respectively; Ra is the terrestrial radiation (MJ m −2 d −1 ).The PET (Eto) units are in mm/year.
1 km.The potential evapotranspiration-PET (Eto)-was calculated as follows: where Tmax and Tmin are annual average maximum and minimum air temperatures (°C), respectively; Ra is the terrestrial radiation (MJ m −2 d −1 ).The PET (Eto) units are in mm/year.

Soil Data and Plant Available Water Content Data
The layer of soil that limits root growth, called the Root Restricting Layer (RRL), is where plant roots stop developing.We used a global dataset map for the depth of this layer from [18], as seen in Figure 4D.
Plant Available Water Content (PAWC) is the amount of water in the soil that plants can consume, calculated by taking the difference between the volumetric field capacity and the permanent wilting point.In our research, we used a global PAWC map from [19].This map shows the available water content for soil up to 200 cm deep, divided into seven layers and combined into one file.As shown in Figure 4C, PAWC values in these river subbasins range from 0.14 mm/mm in the soils of the upper Casanare River subbasin to

Soil Data and Plant Available Water Content Data
The layer of soil that limits root growth, called the Root Restricting Layer (RRL), is where plant roots stop developing.We used a global dataset map for the depth of this layer from [18], as seen in Figure 4D.
Plant Available Water Content (PAWC) is the amount of water in the soil that plants can consume, calculated by taking the difference between the volumetric field capacity and the permanent wilting point.In our research, we used a global PAWC map from [19].This map shows the available water content for soil up to 200 cm deep, divided into seven layers and combined into one file.As shown in Figure 4C, PAWC values in these river subbasins range from 0.14 mm/mm in the soils of the upper Casanare River subbasin to 0.39 mm/mm in the soils of the upper Meta River subbasin, with an average of 0.2 mm/mm in the soils of all studied subbasins.

Land-Use Data, Land-Cover Data, and Kc
The land-use/land-cover (LULC) map used for this study was made by IDEAM, which includes land-use and land-cover data from 2014 to 2018 [22].This map was converted into a raster file simplified into 19 different types of land-use (see Figure 4B).In addition, we created a biophysical table in a comma-separated values (CSV) format.This table consists of five columns and contains detailed information about each LULC category.It includes the land-use (LU) code, a description of the LULC, Kc value, root depth, and the type of vegetation in the LULC categories.Table 2 shows the first three columns of the biophysical table used.The crop coefficient (Kc) is a dimensionless quantity used in this study to calculate crops' water requirements during their growth phases.The FAO (Food and Agriculture Organization of the United Nations) created a well-known set of Kc values for a variety of crops.Ranging from 0 (no water loss) to 1 (maximum water loss), these values are based on studies carried out in different regions and conditions over the world and were obtained from [9].

River Water Discharge Data
This study examines data from five hydrometeorological stations (Table 3).The Aceitico gauging station, located downstream in the study basin, was identified as the main outflow point.The highest annual average water discharge was recorded there at 9288.5 m 3 /s in 2021, while the lowest one at 3647.6 m 3 /s was in 1992.Over the period 1983-2021, the annual average water discharge was 5256.8 m 3 /s at this station.For information on the studied subbasins, see Table 3.This study used climate change data from the WorldClim database, focusing on scenarios available at a 2.5-min resolution [23].These data were downscaled to 1 km resolution using the thin-plate splines interpolation method, adhering to the methodology detailed by [24].The baseline for anomaly calculations was set using the average values from spatial data used for the calibration period from 1983 to 2012.
Key statistical measures-maximum value, minimum value, mean, and standard deviation-for critical climate variables associated with climate models used in this study are presented in Tables 4 and 5.In terms of precipitation and evapotranspiration predictions for the year 2050 under the Shared Socioeconomic Pathways (SSP) Scenario 4.5, the data show considerable variation among the climate models.The BCC-CSM2-MR model predicts the highest maximum precipitation at 5293.1 mm per annum, indicating the likelihood of significant rainfall events.At the lower end, the IPSL-CM6A-LR model records the smallest minimum precipitation at 802 mm per annum, suggesting less precipitation in certain zones of the river basin studied.When considering average conditions, the BCC-CSM2-MR model also leads with the highest mean precipitation at 3115.3 mm per annum.For evapotranspiration, the highest mean is predicted by the HadGEM3-GC31-LL model at 1775.5 mm per annum, while the highest maximum evapotranspiration was also recorded in the HadGEM3-GC31-LL model at 1937.2 mm per annum.
We observe significant variability in precipitation patterns across the models under the Shared Socioeconomic Pathways (SSP) Scenario 8.5.The BCC-CSM2-MR model predicts the highest maximum precipitation at 4904.6 mm per annum, suggesting intense rainfall events could become more common.Conversely, the UKESM1-0-LL model shows the lowest minimum precipitation at 770.3 mm per annum, indicating potential aridity in certain regions.The average precipitation is the highest in the BCC-CSM2-MR model-2837.5 mm per annum, which may imply a general increase in precipitation levels.Regarding evapotranspiration, the UKESM1-0-LL model predicts the highest mean at 1814.8 mm per annum and shows the highest maximum evapotranspiration at 1992.6 mm per annum, reflecting potential water loss from the surface and indicating the possibility of increased drought due to rising air temperatures.

Calibration and Validation
The InVEST-AWY model was calibrated using annual data from 1983 to 2012 obtained from the Aceitico gauging station.The subsequent review phase lasted from 2013 to 2021, as shown in Figure 5.The model's performance was assessed to reduce the average bias and refine the coefficient of determination (R 2 ) along with the Root Mean Square Error (RMSE) (Table 6).Table 6.Metrics for the model performance in the studied subbasins during calibration and validation periods using the InVEST-AWY model (for symbols, see Figure 5).Source: [9].Table 6.Metrics for the model performance in the studied subbasins during calibration and validation periods using the InVEST-AWY model (for symbols, see Figure 5).Source: [9].A total of 100 simulations were carried out annually from 1983 to 2012 [9] with variations in the Kc and Z coefficients.The analysis revealed that the InVEST-AWY model gave optimal results with Z = 1 and Kc = 1.10.Despite identifying the most effective parameter combination for the Meta River basin, only the upper Meta River subbasin showed strong correlation coefficients, scoring 0.79 in the calibration phase and 0.83 in validation, indicating that the model performed accurately in this area.

Results
The simulated average water yield for this basin was 1.62 × 10 11 m 3 /year (5141.6 m 3 /s) (Figure 6A) during the period 1983 to 2021, which is 1.5% higher than the value reported by IDEAM.

Changes in the Annual Water Yield Predicted 2050 under CMIP6 Scenarios
Table 7 presents annual water yield data (in m 3 /year) for various subbasins and scenarios, including the current baseline (1983-2012) and predictions for the year 2050 under two Shared Socioeconomic Pathways (SSP 4.5 and SSP 8.5 scenarios).In addition, the percentage change in annual water yield volume and flow changes for these future scenarios compared with the baseline period was revealed.

Changes in the Annual Water Yield Predicted 2050 under CMIP6 Scenarios
Table 7 presents annual water yield data (in m 3 /year) for various subbasins and scenarios, including the current baseline (1983-2012) and predictions for the year 2050 under two Shared Socioeconomic Pathways (SSP 4.5 and SSP 8.5 scenarios).In addition, the percentage change in annual water yield volume and flow changes for these future scenarios compared with the baseline period was revealed.The Meta River basin shows a slight increase in annual precipitation under both future scenarios (Table 7), with a corresponding rise in PET.The simulated water flow increases according to the SSP 4.5 scenario but reveals a smaller rise according to the higher emission SSP 8.5 scenario (Figure 6; Table 7).Water yield changes reflect an increase of 24% under the SSP 4.5 scenario and 19% under the SSP 8.5 scenario, indicating potential water resource availability.
The upper Meta River subbasin experiences minimum changes in annual precipitation across the scenarios (Table 7), while PET shows a moderate increase.The simulated water flow and water yield changes under future scenarios indicate a slight decrease under the SSP 4.5 scenario with −1% and under the SSP 8.5 scenario with −6% (Figure 6; Table 7).In the case of the upper Casanare River subbasin, precipitation and PET are predicted to increase under both scenarios.Remarkably, the changes in water yield are noticeable, with a significant increase of 15% under the SSP 4.5 scenario and 9% under the SSP 8.5 scenario (Figure 6; Table 7).It is also important to note the InVEST-AWY overestimation in 108% in this subbasin, which means careful recalibration for that specific zone is necessary.
For the Yucao River subbasin, the model predicts an increase in annual precipitation and a slight increase in PET for both future scenarios (Table 7).The changes in water yield are markedly positive, with increases of 18% and 14% under the SSP 4.5 and SSP 8.5 scenarios (Figure 6; Table 7), respectively, indicating a trend towards greater water resource availability.
Predictions for the South Cravo River subbasin show decreased annual precipitation and minor changes in PET under both future scenarios.The simulated water flow will decrease significantly, with a decrease of −10% under the SSP 4.5 scenario and −14% under the SSP 8.5 scenario.This prediction suggests a concerning decline in water availability for the South Cravo River subbasin.However, the model performance for the South Cravo River subbasin is notably less accurate with a sub estimation of 41% in the baseline period, indicating a need for model refinement.The predicted negative impact and the discrepancy with observed data emphasize the need for improved data and modeling at the local level.
Based on the InVEST-AWY model's predictions for the year 2050 under the SSP 4.5 and 8.5 scenarios, there is an anticipated increase in available water in 24% and 19%, respectively, throughout the entire basin of the Meta River, attributed to heightened rainfall in the higher elevations.Simultaneously, scenarios [25] suggest a substantial decrease, up to 85%, in low water flows in over half of the Colombian Orinoco River basin.
Table 7 also highlights the variability and uncertainties inherent in the InVEST Annual Water Yield model's predictions when calibrated with data from 1983-2012.The discrepancies in water yield changes are considerable, with the model overestimating yields by 108% in the Upper Casanare River subbasin, indicating a possible overprediction of water availability in this area.Conversely, the South Cravo River subbasin's yield is underestimated by 41%, suggesting an underprediction of water yield.These errors could be attributed to the sparse network of weather stations within these watersheds and the significant topographical variations, both of which are challenging factors for the model to accurately simulate [16].

Spatial Variation in Water Yield for 1983-2012 and 2050 under Two Scenarios
Figure 6 shows the outputs from the annual water yield model identifying the fluctuating availability of water across a historical baseline (annual average water yield for the period 1983-2012) and two prediction scenarios for 2050, based on CMIP6 SSP 4.5 and SSP 8.5 scenarios.The historical baseline, covering the period from 1983 to 2012, frames the analysis, with a noticeable gradient in water yield across the landscape.Higher elevation zones located in the western part of the Meta River basin, typically associated with higher annual precipitation, present a reduction in water yield.This decline is attributed to fluctuations in precipitation coupled with rising air temperatures, which are evident in the shift from historical patterns to the 2050 predictions.The SSP 4.5 scenario, considered moderate, retains much of the baseline's spatial pattern, but it shows a subtle shift in the distribution of water yield across these higher elevations, signaling a nuanced response to evolving climatic conditions.
On other hand, the SSP 8.5 scenario, characterized by higher greenhouse gas emissions, predicts a severe deviation from historical trends, especially in the western highlands, where more intense changes are expected.This scenario indicates an increased change in water yield, and a deeper color on the map (Figure 6) signals a potential increase in water availability.Such changes are likely a consequence of predicted increases in precipitation and the amplifying effects of increased air temperatures on evapotranspiration processes.Meanwhile, low-lying areas of the Meta River basin are predicted to experience increased rainfall.

Model Limitations
The InVEST-AWY model uses annual average data and does not consider the changes in water supply and hydropower production that occur at different times of the year, such as dry or rainy seasons.This means that predicting the amount of water or energy available during extreme weather conditions such as drought or flooding may not be very accurate.Also, the way the model calculates how much water is used up (consumptive demand) is too simple, which may affect how well it can determine how much water is available for various purposes.This is important when looking at how land-use and land-cover changes affect water resources [17].Furthermore, the InVEST-AWY model does not separate surface and groundwater runoff [26].On the other hand, the model shows a minimum response to variations in the Z coefficient [27,28], but it is very sensitive to changes in precipitation and evapotranspiration [29].

Uncertainties in the Reported Information
The lack of comprehensive hydrogeological studies in the Meta River basin makes it difficult to fully understand hydrogeological dynamics, especially the contribution of groundwater to total river flow [29].This finding points to diverse water dynamics within aquifers, suggesting more complex interactions between surface and groundwater systems than previously thought.
On the other hand, this study faces significant uncertainty due to incomplete data from the selected monitoring stations.Among the five stations used, only one is automated, and notable data gaps were observed, such as missing annual data for certain years at the Aceitico and Cravo Norte gauging stations, as well as Puente Yopal (see Figure 5).In addition, these stations reported irregular values in various years, potentially impacting the reliability of the model calibration and validation.This circumstance highlights the difficulty of accurately estimating water flow in this area.Furthermore, Pimentel et al. [25] found that although groundwater extraction accounts for approximately 5% of total regional demand, the actual figure may be higher due to unaccounted-for extraction.

Climate Change Uncertainties
In the context of the 2050 water yield predictions under the CMIP6 for the SSP 4.5 scenario, Figure 7 shows clear spatial differences in predictions between different global climate models.The BCC-CSM2-MR model predicts a significantly high-water yield, especially in the western part of the study basin, with estimates ranging from 3000 to 4500 mm per annum.This increase is particularly noticeable on the mountain slopes in the western region, in the southwest of the basin.In contrast, models such as HadGEM3-GC31-LL, MPI-ESM1-2-HR, and UKESM1-0-LL indicate a drier trend in most parts of the basin.However, on the mountain slopes of the river basin, they do predict isolated zones with maximum water yields of about 1500 mm per annum.The remaining models generally suggest a moderate increase in water yield in the mountainous areas, maintaining a trend that falls between these two extremes.
Hydrology 2024, 11, x FOR PEER REVIEW 14 of 1 especially in the western part of the study basin, with estimates ranging from 3000 to 450 mm per annum.This increase is particularly noticeable on the mountain slopes in th western region, in the southwest of the basin.In contrast, models such as HadGEM3 GC31-LL, MPI-ESM1-2-HR, and UKESM1-0-LL indicate a drier trend in most parts of th basin.However, on the mountain slopes of the river basin, they do predict isolated zone with maximum water yields of about 1500 mm per annum.The remaining models gener ally suggest a moderate increase in water yield in the mountainous areas, maintaining trend that falls between these two extremes.Regarding the uncertainties inherent in these climate predictions, the BCC-CSM2 MR model exhibits the highest level of variance.This is primarily due to its significan deviation from the ensemble average (estimated from all 13 models used in this study with an overestimation ranging from 400 to 1200 mm, especially noticeable on the moun Regarding the uncertainties inherent in these climate predictions, the BCC-CSM2-MR model exhibits the highest level of variance.This is primarily due to its significant deviation from the ensemble average (estimated from all 13 models used in this study) with an overestimation ranging from 400 to 1200 mm, especially noticeable on the mountain slopes of the study basin.At the other end of the spectrum, models like FIO-ESM-2-0, HadGEM3-GC31-LL, and IPS-CM6A-LR show strong underestimations, deviating from the ensemble average by approximately −150 to −450 mm.Notably, the MIROC6 model demonstrates the least uncertainty, maintaining water flow deviation in a narrow range from 0 to 100 mm.
In the high-emission SSP 8.5 scenario (depicted in Figure 8), the CMIP6 model water yield predictions show a pattern largely similar to the SSP 4.5 scenario, with some notable variations.An analysis of the map mosaic shows that most models predict a small increase in water yield, mainly concentrated on the mountain slopes of the study basin.This trend is consistent across various models, with the BCC-CSM2-MR model continuing to predict high water yields.In contrast, models such as HadGEM3-GC31-LL, MPI-ESM1-2-HR, and UKESM1-0-LL predict lower to moderate water yields.
Hydrology 2024, 11, x FOR PEER REVIEW 1 variations.An analysis of the map mosaic shows that most models predict a small in in water yield, mainly concentrated on the mountain slopes of the study basin.This is consistent across various models, with the BCC-CSM2-MR model continuing to p high water yields.In contrast, models such as HadGEM3-GC31-LL, MPI-ESM1-2-H UKESM1-0-LL predict lower to moderate water yields.In terms of uncertainties, the BCC-CSM2-MR model again demonstrates a sign overestimation of annual water yield under the SSP 8.5 scenario, 150 to 900 mm abo average of all the models.This overestimation underscores the model's tendency dict higher water yields, especially under higher emission scenarios.In addition, t Earth3-Veg model predicts a substantial increase in water yield, particularly in the c western zone of the river basin, with estimates ranging from 250 to 500 mm.The ESM2-0 model also indicates a modest positive deviation in water yield for the nort ern zone in the range of 250 to 300 mm.In contrast, models such as INM-CM5-0 an In terms of uncertainties, the BCC-CSM2-MR model again demonstrates a significant overestimation of annual water yield under the SSP 8.5 scenario, 150 to 900 mm above the average of all the models.This overestimation underscores the model's tendency to predict higher water yields, especially under higher emission scenarios.In addition, the EC-Earth3-Veg model predicts a substantial increase in water yield, particularly in the central-western zone of the river basin, with estimates ranging from 250 to 500 mm.The MRI-ESM2-0 model also indicates a modest positive deviation in water yield for the north-eastern zone in the range of 250 to 300 mm.In contrast, models such as INM-CM5-0 and MIROC6 show the least uncertainty, with deviations ranging from 0 to 150 mm across most of the river basin.
In the most severe of these scenarios, the Meta, Vichada, and Guaviare rivers are expected to experience significant reductions in low water flow, by 95%, 98%, and 50%, respectively.Moreover, in a study carried out by [16], using the SWAT model, simulations were carried out in the Meta River basin, covering the basins of the Ariporo, Cravo Sur, Cusiana, Guatiquía, and Guayuriba rivers.This study aimed to assess the impact of various climate change scenarios (RCP 2.6, RCP 4.5, and RCP 8.5 scenarios) for the year 2040, using a reference period of 1997 to 2012 (16 years).These predictions modeled annual water flow reductions of 6-7%, with January identified as the most critical month for water flow reductions.

Conclusions
(1) For the Meta River basin, we predict a significant increase in simulated water flow from the current 5141.6 m 3 /s to 6397.5 m 3 /s by 2050 under the SSP 4.5 scenario, and an increase to 6101.5 m 3 /s under the SSP 8.5 scenario.This correlates with an increase in water yield by 24% and 19%, respectively, under the two future scenarios evaluated.The upper Meta River subbasin shows a slight decrease in water flow from the current 1600.6 m 3 /s to 1578.3 m 3 /s (SSP 4.5) and a decrease to 1511.3 m 3 /s (SSP 8.5), with water yield changes ranging by −1% and −6%, respectively.The upper Casanare River subbasin is expected to see a moderate rise in water yield from 976.3 m 3 /s to 1119.7 m 3 /s and 1059.9 m 3 /s under the SSP 4.5 and SSP 8.5 scenarios, respectively, with water yield changes rising by 19% and 9%, respectively.The Yucao River subbasin shows an increase from 95.9 m 3 /s to 113.3 m 3 /s and 109 m 3 /s, with water yield changes increasing by 18% and 15%, respectively.In contrast, the South Cravo River subbasin is predicted to face a decrease in water flow from the current 58.5 m 3 /s to 52.9 m 3 /s and 50.3 m 3 /s, with a significant drop in water yield changes of −10% and −14%, respectively, indicating a marked reduction in water availability.(2) Although the InVEST-AWY model provided acceptable results for the entire Meta River basin using data from 1983 to 2012, our study showed that the model is capable of effectively predicting potential impacts in well-calibrated areas, especially in the upper Meta River subbasin as defined by the Humapo gauging station.(3) The uncertainties observed in the thirteen global climate models according to SSP 4.5 and 8.5 scenarios stem from the varying predictions of increased water yield availability in the flatter regions of the main basin.This potential increase could lead to a higher likelihood of concurrent floods or river overflows, emphasizing the need for adaptation strategies in these areas.(4) Future research should prioritize two key areas.First, flood risk analysis and strategies are needed in areas that have potential for increased water yield, considering expected increases in water levels and the possibility of flooding.This area of research will include the use of models to predict floods, assess impacts on infrastructure and communities, and develop flood mitigation strategies.Second, it is important to study the socioeconomic impacts of water yield fluctuations, especially in regions facing declining water availability, such as the South Cravo River subbasin.Research in this area may focus on water management, impacts on agricultural practices, and impacts on community livelihoods.(5) Finally, it would be useful to carry out comparative studies using a range of both nonrobust and robust hydrological models.This approach would serve to validate the study's findings and provide a clearer understanding of the comparative advantages

Figure 1 .
Figure 1.The location of the Meta River basin (Research zone) and its subbasins in Colombia.

Figure 2 .
Figure 2. Flowchart for calculating changes in water yield in the Meta River basin using the InVEST-AWY model.

Figure 1 .
Figure 1.The location of the Meta River basin (Research zone) and its subbasins in Colombia.

Figure 1 .
Figure 1.The location of the Meta River basin (Research zone) and its subbasins in Colombia.

Figure 2 .
Figure 2. Flowchart for calculating changes in water yield in the Meta River basin using the InVEST-AWY model.

Figure 2 .
Figure 2. Flowchart for calculating changes in water yield in the Meta River basin using the InVEST-AWY model.

Hydrology 2024 , 18 Figure 5 .
Figure 5. Fluctuations in observed (represented by a solid blue line) and modeled (depicted with a solid red line) water discharge (Q) in the studied rivers of the Meta River basin from 1980 to 2021.The dotted line traces the sixth-degree polynomial trend of water discharge, whether observed or modeled.R 2 denotes the fit quality of the trend line to the data points; rcal and rval signify the correlation coefficients for the calibration (1983-2012) and validation (2013-2021) period, respectively.Source: [9].

Figure 5 .
Figure 5. Fluctuations in observed (represented by a solid blue line) and modeled (depicted with a solid red line) water discharge (Q) in the studied rivers of the Meta River basin from 1980 to 2021.The dotted line traces the sixth-degree polynomial trend of water discharge, whether observed or modeled.R 2 denotes the fit quality of the trend line to the data points; r cal and r val signify the correlation coefficients for the calibration (1983-2012) and validation (2013-2021) period, respectively.Source: [9].

Figure 7 .
Figure 7. Simulated annual water yield in the Meta River basin and associated uncertainties for eac global climate model according to the scenario SSP 4.5.

Figure 7 .
Figure 7. Simulated annual water yield in the Meta River basin and associated uncertainties for each global climate model according to the scenario SSP 4.5.

Figure 8 .
Figure 8. Simulated annual water yield in the Meta River basin and associated uncertainties f global climate model according to the SSP 8.5 scenario.

Figure 8 .
Figure 8. Simulated annual water yield in the Meta River basin and associated uncertainties for each global climate model according to the SSP 8.5 scenario.

Table 2 .
Crop coefficient (Kc) used for each LULC category.

Table 3 .
The gauging stations used in the study (AAWD-annual average water discharge).

Table 4 .
Descriptive statistics for the models according to the SSP 4.5 scenario for the Meta River basin.

Table 7 .
Spatial variation in the annual water yield in the Meta River basin under the SSP 4.5 and 8.5 scenarios.
Note: PET is Potential Evapotranspiration, AET is Actual Evapotranspiration; 1 uncertainty estimated in the baseline simulation for the period 1983-2012.