Hydrological Modeling for Multifunctional Landscape Planning in the Orinoquia Region of Colombia

With over 200,000 km2 of natural savannas, the Orinoquia region of Colombia is a key and strategic conservation area. Because of Colombia’s fast economic growth, there are significant plans for agro-industrial expansion in the Orinoquia. This expansion may seriously affect water availability. To evaluate the cumulative impacts on freshwater ecosystems derived by different expansion scenarios, the use of a comprehensive framework for mathematical modeling, able to represent the hydrological processes at a macro-basin scale, is crucial for analysis and as a tool to bridge the gap between science and practice. In this work, we developed a general methodological framework for hydrological analysis at macro-basin scale consisting of four main stages: 1) collection and processing of hydro-climatological data, 2) characterization of hydro dependent water use sectors, 3) mathematical modeling and 4) scenario simulation. As a result of applying the proposed framework, we obtained a coupled hydrological model, which allows us to represent the rain-runoff process, the river-floodplain interaction and anthropic processes such as surface water extraction and groundwater extraction, enabling us to represent the complexity of the Orinoquia region. The model was successfully implemented in Matlab showing a Nash-Sutcliffe efficiency coefficient between 0.62 and 0.92 in calibration and between 0.49 and 0.92 in validation. With this model we analyzed five different agro-industrial expansion scenarios, finding that the Colombian Orinoquia may have future high pressure on water resource areas with critical changes in the water availability regime. The scenarios show reductions of up to 85% in low water flows in more than 50% of the area of the Colombian Orinoco basin. In the most extreme scenarios, the Meta, Vichada and Guaviare rivers show reductions of 95, 98 and 50% in low water flows. The results show an urgent need to consider hydrology in planning processes to ensure the sustainability of this important area in Colombia.


INTRODUCTION
Traditionally, extensive cattle ranching has played an important role in the economy of the Colombian Orinoco (Andrade and Castro, 2012;Buriticá, 2016). However, in the last 2 decades, agriculture has expanded, reaching a significant increase from 1,000 km 2 to 8,000 km 2 of planted area (around 2.3% of the Colombian Orinoco basin). This expansion of the agricultural Frontier (Vargas et al., 2015;Cárdenas-Santamaría et al., 2018) may have a strong impact in savanna ecosystems, derived from the alteration of the land for crops such as rice, sugar cane and oil palm (Andrade and Castro, 2012). At the same time, it has led to an increase in water consumption for the different crops, generating substantial pressure on the water resources of this region.
The National Planning Department of Colombia (DNP, 2018) estimated that more than 150,000 km 2 of the Colombian Orinoquia is suitable for agriculture, suggesting greater pressure on the water resources of this region in the future. Globally, agriculture consumes 70% of the global freshwater withdrawals due to crop evapotranspiration (FAO, 2017). Similarly, a large increase in crop area in the Orinoquia region would certainly affect water quantity and quality. Research has clearly established how agricultural expansion can have severe consequences for ecosystems (Byard, 1999;Kanianska, 2016;Rabin et al., 2020;Tamaris et al., 2017). The Colombian Orinoquia covers 8.5% of the flooded savannas in South America and it is a strategic ecosystem of great economic, biological and ecological importance for the country (Mora et al., 2015). It is also one of the most important biodiversity reserves in the neotropics (Gassón, 2002).
Adequate planning of the territory must be based on multiscale analysis, supported by measurements and mathematical modeling tools. These tools allow for prospective analysis and evaluation of the impacts on the territory in the face of future interventions. In Colombia, historically, most of the hydrological modeling efforts have focused on the Magdalena-Cauca basin, considered the most important basin in the country (Angarita et al., 2018;Craven et al., 2017;Labrador Cadena et al., 2016;Rodríguez et al., 2020). However, the Orinoquia region has been identified as the next major area of agro-industrial expansion in Colombia, so it is urgent to develop tools for prospective analysis. Considering this need, this study presents a hydrological modeling exercise with monthly time step for the Colombian Orinoquia region, addressed by two parsimonious conceptual models: the ABCD model proposed by Thomas (1981) and the river-floodplain interaction model proposed by Angarita et al. (2018). We developed a logical technical framework allowing for the description and understanding of the territory that, in a practical way, enables the evaluation of possible future scenarios and their impacts on water availability.

STUDY AREA
The Colombian Orinoquia region ( Figure 1) has an area of 347,208 km 2 covering the departments of Arauca, Casanare, Vichada, Guainía, Cundinamarca, Boyacá, Norte de Santander, Vaupés, Meta, Guaviare and Santander; it represents about 34% of the Colombian territory. It exhibits a unimodal pattern in the temporal distribution of precipitation in the annual cycle, with an average of 2,650 mm yr −1 , which is largely associated with the greater convective activity generated by the southern migration of the Intertropical Convergence Zone (Poveda et al., 1998;Poveda and Mesa, 1996;Poveda and Rojas, 1997;Poveda et al., 2002;J. Velez et al., 2000). Temperature can vary between 23°C and 29°C during the year. Geologically, it consists of tertiary sedimentary rocks, which cover crystalline Precambrian rocks and sediments of the Paleozoic and Cretaceous (SGC, 2020). It also has a large variety of geomorphological structures that allow for the presence of a very diverse and complex drainage system, ranging from meandering to braided rivers. Thanks to the low slopes of its plain areas, the broad transversal sections of the channels (around 0.5-2 km) and the constant presence of land depressions around the rivers, there are bidirectional interactions between rivers and their savanna areas, generating large floods.
Consistent with the rainfall regime, river flows are highest in June, July and August. The main tributaries that are part of the Orinoco basin in Colombia are the Meta river with a multi-year annual average flow of 4,973 m 3 s −1 , the Guaviare river with 4,035 m 3 s −1 , the Inírida river with 2,978 m 3 s −1 and the Vichada river with 1,030 m 3 s −1 . In total, the Orinoco basin in Colombia provides 26.4% of the fresh water of the Colombian territory (IDEAM, 2019).

METHODOLOGY
The study of water systems implies knowing their character, behavior, distribution and condition (Brierley and Fryirs, 2004). Therefore, for this study we have conceptualized four main stages to develop the analysis: 1) collection and processing of hydroclimatological data, 2) characterization of hydro dependent water use sectors, 3) mathematical modeling (this step includes calibration and validation of the model) and 4) scenario simulation. Each of the stages is described below.

Collection and Processing of Hydro-Climatological Data
The Colombian territory has a large hydro-climatological monitoring network (over 2,600 active stations). However, most of the stations are in the Magdalena-Cauca basin. While the Orinoco basin represents 34 percent of the national territory, it contains only about 14% of all active monitoring stations in Colombia. Most of the stations are located on the eastern mountain range, leaving the savanna areas (more than 70%) poorly monitored. We compiled the hydrometeorological information used in this study from the stations of the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM). In total, 526 monitoring stations were identified for both meteorological and hydrological variables. Of the 526 stations identified in the study area, information was acquired from 468 of them, with a valid registration period until 2014. Most of the stations were on the western side of the study area, while the central and eastern regions presented low station coverage. For all time-series data, we identified anomalous data, analyzed consistency and homogeneity, and replaced missing data. Processes carried out are described below.
Identification of Anomalous Values: We detected anomalous data using three methodologies: the Grubbs test, Tukey's test, and the test of double absolute deviations from the median. For each variable, we analyzed the detection rate and correlation with climatic variability. The Oceanic Niño Index (ONI) correlated the most with each variable.
Consistency analysis: After the anomalous data analysis, our approach tested completeness, considering the expected and observed length of the series. As a first criterion, the effective period of registration of the series cannot be less than 10 years. Additionally, we included a consistency criterion, defined as less than 20% of missing data. To pass the consistency analysis, the series had to meet both criteria simultaneously.

Estimation of missing data
Estimation of missing data was carried out by combining the Inverse distance weighting (IDW) and Ordinary Kriging methodologies for precipitation and multiple linear regressions for temperature.

Homogeneity
The mass curve method was applied for all variables of interest. Thus, the station with the highest number of records, the least amount of missing data and the least estimate of anomalous data was taken as the reference station. If a series is homogeneous, it will present linear behavior in relation to the reference station. If not, it will present significant changes in the slope.
According to the analysis of consistency and homogeneity of the information, 227 precipitation stations, 59 temperature stations and 57 flow stations were considered acceptable.

Characterization of Hydro Dependent Water Use Sectors
For our analysis, we considered five hydro dependent water use sectors: 1) agriculture, 2) livestock, 3) domestic, 4) hydrocarbons, and 5) mining.
We established the use of water in agricultural production based on the irrigation needs of the different crops by evaluating the amount of water and the time of its use. This approach achieved a balance between the amount of water required by the crop, in compensation for the loss by evapotranspiration, and the effective precipitation. In this sense, the water demand for the agricultural sector is considered as the water extracted to supply the water requirements of the crops through irrigation. The Colombian Ministry of Agriculture (MinAgricultura, 2018) provided the information for the estimate from the Municipal Agricultural Assessments for the years 2014 and 2015.
The demand of the livestock sector was based on the water required to satisfy the basic needs of the animal during its life cycle, as well as the water consumption required in the slaughter and benefit stages for cattle, pigs and birds. The base information used to estimate demand in this sector corresponded to 2016 National Livestock Census (ICA, 2016).
For the domestic sector, we considered the water demand as the water used in household activities. These may include direct drinking, food preparation, personal hygiene and cleaning of elements, materials, or utensils. The population information taken for the calculations was the projections made by DANE in the 2005 National General Census, with the base year 2016 (DANE, 2008).
For the hydrocarbons sector, we estimated the demand considering all stages of production, transportation and refining of crude oil. The base information considered for estimating demand in this sector was the production of crude oil during 2016 consolidated in the Petroleum Statistical Report (ACP, 2017).
For the mining sector, we estimated the amount of water required to produce 1 cubic meter or ton of mineral. The base information considered for estimating the demands of the mining Frontiers in Environmental Science | www.frontiersin.org October 2021 | Volume 9 | Article 673215 sector was the mineral production report of the National Mining Agency (ANM, 2018) for 2015. All surface water demands were estimated using the methodology described by Bedoya, Contreras, and Ruiz (2010). The surface water returns were considered as a fraction of the surface water demand as well as the ground water demand. For this case, a return of 80% was considered for the domestic sector and 50% for the livestock sector. We did not include returns in the other sectors.
Groundwater extraction supplies about 5% of the total demand from the sectors, according to the reports of licenses granted by the environmental authorities. However, it is possible that this value is higher, due to illegal extractions or those that are otherwise unreported to authorities. Given the uncertainty that exists in monitoring groundwater demands, this factor was considered in the model as an effective calibration parameter.

Conceptual Model
Most of the hydrological phenomena and sub-processes involved in the rainfall-runoff process take place in sub-daily resolutions, so their representation requires the use of models designed in similar resolutions. However, for this study, we conceptualized a tool for the global estimation of the hydrological dynamics of the Colombian Orinoquia, in such a way that a monthly analysis was adequate.
Because the Colombian Orinoco basin is poorly monitored, we built a model requiring little information for its construction and operation, with parsimonious parameterization in its equations. For the conceptual model, it is important to consider that the Orinoco basin has about 480,000 km 2 of natural savannas, representing 17.8% of the natural savannas in South America (Mora et al., 2015). Of these, 12.5% correspond to flooded savannas. The flood pulse is the main driving force that modulates the annual changes in biotic and abiotic variables that take place in the main channel and in all water bodies associated with the floodplain (Junk et al., 1989). The high geomorphological complexity of the flooded savannas results in high heterogeneity of habitats and high biological productivity; broad biodiversity is maintained over time thanks to the action of periodic floods (Mora et al., 2015). Therefore, our approach includes a representation of the interactions between the river and the natural savannas.
We also included the representation of surface water and groundwater demands in the conceptualization of the model. Inclusion of groundwater is important because of the National Government´s notable interest in expanding the hydrocarbons, agriculture, livestock, mining and domestic sectors in the Orinoco basin (DNP, 2018), Our technical analysis of the characteristics of the Colombian Orinoquia suggested the following main hydrological processes for the monthly time step representation: Evaporation, transpiration, precipitation, base flow, infiltration, deep percolation, surface water demands, groundwater demand, and the bidirectional interactions between rivers and natural savannas. Figure 2 shows the scheme of the conceptual model proposed for the analysis in this study of the Orinoco basin in Colombia.

Mathematical Model and Computational Algorithm
In accordance with the conceptual model, we reviewed 14 well known hydrological models to determine the most appropriate mathematical model in line with the scope of the research (See supplementary material). As a result of this review, we considered the ABDC hydrological model coupled with the floodplain model proposed by Angarita et al. (2018) as the best option that allows us to represent most of the conceptualized processes through a simple and parsimonious scheme. The representation of the river floodplain interaction phenomenon is highly complex. Studies such as the one carried out by Pontes et al. (2017) show a representation of the complexity of this process based on an October 2021 | Volume 9 | Article 673215 inertial model. In the case of the present study, the simplified conceptual model presented by Angarita et al. (2018) was chosen in order to have a simple approach to the problem. The coupled model was configured under a semi-distributed spatial calculation scheme, considering the Hydrological Analysis Unit (HAU) as the minimum unit of analysis. To couple the ABCD model with the Angarita et al. (2018) river-floodplain interaction model, a cumulative calculation scheme was used. This implies that, at each step of time i, the water balance is solved in the j Hydrological Analysis Unit (HAU j ) using the equations of the Thomas (1981) model.
The ABCD model (Liu, 2020;Wang et al., 2020;Zhang et al., 2020) considers two water storage components for the analysis of the water balance in the HAU j ; the first component represents the soil or evapotranspiration zone (Sw i,j ) above the water table and the second component represents the saturated zone ((Sg i,j ) below the water table.
For modeling purposes, the subsurface flow component in the surface part of the evapotranspiration zone can be included in the direct runoff (Ro i,j ). The model considers negligible the lateral deep flow in the unsaturated zone, in such way that the potential recharge is equaled to the real recharge (Rg i,j ).
Applying the continuity equation to a control volume Sw i,j : Where: P ij , precipitation; ET i,j , actual evapotranspiration; Rg i,j , recharge; Ro i,j , direct runoff; Sw i,j , change in soil storage between the calculation period i and the immediately preceding period (Sw i−1,j ). Thomas (1981) further defines the variables W i,j (available water) and Y i,j as: In each time interval, we assume that the humidity decreases according to the exponential decay law, assuming as initial humidity at the beginning of each interval Y i,j , where ETP i,j is the potential evapotranspiration and b j is a parameter of the model: Thomas (1981) defines the state variable Y i,j as a non-linear function of water available according to the parameters a j (dimensionless) and b j : Where a j and b j are parameters that can be determined from measurements of precipitation, evapotranspiration, and soil moisture in the basin. The upper limit of Y i,j is represented by the parameter b j . Thomas (1981) notes that the function Y i,j has no particular meaning. Then, when replacing in the previous equations, we obtain: To differentiate the runoff from the recharge, a partition coefficient c j is assumed: The underground flow Qg i,j , which is the fraction of the flow observed in the river that comes from the storage in the saturated zone (Sg i,j ), is: The storage Sg i,j is interpreted as a dynamic storage that represents the connectivity between the river and the aquifer. Therefore, when applying the continuity equation to a volume of aquifer storage Sg i,j : Where ΔSg i,j is the change in storage of the saturated zone and Sg i−1,j is the storage of groundwater in the immediately preceding period. By substituting the previous equations and solving for Sg i,j : The total flow, that is the flow observed in the river, is: According to Thomas (1981), the parameters a, b, c, and d can be interpreted in the following way: a is the tendency of runoff that occurs before the soil is completely saturated. Values of " a less than one generate runoff when W i < b.
b is the upper limit to the sum of the real evapotranspiration and soil moisture.
c is the fraction of the average flow of the river that comes from the groundwater. d is the reciprocal of the residence time of the groundwater.
Having solved all the previous equations for the time step i, the flows, surface water demands, and water returns of each HAU j are accumulated. As the model accumulates, it subtracts the surface water demands from the river flow and adds the water returns. Therefore, the total flow of the river (Q t i,j ) would be given by: Where Dsp i is the surface water demand and Rsp i the surface water return.
Then, the model checks if there are zones with riverfloodplains interaction. If there is, the model solves the equations of the Angarita et al., 2018 model as follows: Where Q T i,j corresponds to the accumulated flow upstream, Q l i,j represents the lateral flow between the river and the floodplain, R l i,j represents the lateral flows between the river and the floodplain, A j in the area of the HAU j , A fp j is the floodplain area of HAU j , V fp i,j is the volume of water that the floodplain can store, Q thresholdj and V thresholdj indicate the percentage of contribution in each of the directions of the interaction, T rfj and T frj indicates the percentage of contribution in each of the directions of the river-floodplains and floodplains-river interaction. The model leaves groundwater demand unresolved at the HAU level, requiring the use of Hydrogeological Units of Analysis (UHG). A UHG was defined as the grouping of HAUs belonging to the same hydrogeological unit. The reason groundwater demands were managed this way is because most of the time there is no correspondence between the limits of the surface basins and the limits of the aquifers.
For the water balance in groundwater storage, incorporating groundwater demands, the volume of water stored in the UHG (Vsb i ) is calculated as the sum of the groundwater storages of the n HAU that are part of the same UHG.
Where k denotes that a HAU is part of the same UHG and Ext is the percentage of water demand that is supplied by groundwater.

Model Configuration for the Colombian Orinoquia
The modeling framework used the Arc Hydro Tools package in ArcGIS10.3, generating 4,273 HAUs with an accumulation threshold of 10 km 2 . Computationally, the model interprets a HAU as a modeling node and represents the basin segment as a graph. Thus, the framework assigned unique codes to each HAU, establishing connectivity between neighboring HAUs. We assigned the climatic and water demand variables to each unit individually. In the period from 1985 to 2014, the precipitation and temperature fields were generated at monthly temporal resolution, with a pixel size of 250 m. We used Ordinary Kriging methodology for precipitation and multiple linear regressions for temperature to supply this data. In total, 360 raster files were generated. As a synthesis, Figure 3 shows the averages of each variable.
The potential evapotranspiration was estimated using the Thornthwaite and Mather (1957) methodology. The water demands were integrated individually for each sector.
For the evaluation of the scenarios, we assumed that the climatic conditions remain constant in time and space, evaluated with the climatic series of the baseline condition (historical record).

Calibration and Validation
We used the Shuffled Complex Evolution (SCE-UA) heuristics proposed by Duan et al. (1994) to calibrate the hydrological model. Thus, a total of 57 flow measurements were used whose series were verified to be homogeneous and consistent (Figure 4).
In total, 9 parameters (a, b, c, d, Ext, Q threshold , V threshold , T rf , T fr ) were calibrated in each HAU with the presence on floodplains and 5 parameters (a, b, c, d, Ext) in the HAU remaining. For the segmentation of the calibration and validation periods, the traditional approach of 70% of the data to calibrate and 30% to validate was used. However, in the stations that did not have sufficiently extensive series, no validation was performed, so the entire data was used for calibration. The time periods used for each phase can be consulted in the supplementary material. The performance metric that was used in the calibration process was the Nash-Sutcliffe coefficient (NSE). However, metrics such as Absolute Mean Error (AME), Differential of Peak Events (PDIFF), Mean Absolute Error (MAE), Mean Square Error (MSE), Root Mean Square Error (RMSE), Fourth Root of the Mean Error Raised to the Fourth Power (R4MS4E), Root of Absolute Error (RAE) and Percentage Error in Peaks (PEP) were estimated (Teegavarapu and Elshorbagy, 2005;Dawson et al., 2007). Table 1 shows the performance metrics evaluated.
Given the nature of the hydrological model, the calibration was carried out under a cumulative scheme, homogeneously calibrating the HAUs up to the closure point with a flow gauge. Figure 5 presents a conceptual example of the process used. It shows the HAUs associated with three flow gauges.
Knowing that flow gauges one and two are of order 1 (they do not have upstream stations) and flow gauge three is of order 2 (it has one upstream station in line on the same main channel), the calibration process starts hierarchically from the HAUs associated with gauges of order 1, to those of order n. Note that in the example presented in Figure 5, HAU-1, HAU-2 and HAU-3 share the same set of parameters, since they have flow gauge 2 as the gauge in the closure point that can bring flow measurement data for the calibration of these three HAUs. For the proposed model, there were 12 hierarchy orders associated with the available flow gauges. In the HAUs where there were no flow stations, the parameters were regionalized by similarity of morphometric characteristics with calibrated HAUs. Figure 4 shows the HAU groups related to the flow stations used in the calibration. Uncalibrated HAUs are marked in gray.
As part of the analysis, we proposed an evaluation of the model exploring its performance, the suitability of its structure, the identifiability of its parameters, and the uncertainty of its predictions. To this end, 10,000 Monte-Carlos simulations were carried out using the Latin hypercube sample (LHS) method and the responses of the model were analyzed with the Monte Carlo Analysis Toolbox (MCAT) (Wagener et al., 2002) using dotty plots, regional sensitivity analysis and Generalized Likelihood Uncertainty Estimation (GLUE) (Beven and Binley, 1992).  p100

Percentage Error in Peaks
In the following equations, Qi is the observed value,Qi is the modelled value (where i 1 to n data points) and Q is the mean of the observed data Frontiers in Environmental Science | www.frontiersin.org October 2021 | Volume 9 | Article 673215 8 historical commodity trends, existing policies, and stakeholder views, ranging from least to greatest intervention. Scenario 1 may be considered as business as usual, representing a modest expansion that is optimally generated in terms of maximizing agricultural benefits. Scenario two represents the vision of small farmers, ranchers, and agricultural industry on Colombia Orinoquia region. Scenario three represents the vision of the Colombian government described in the Master Plan (DNP, 2016), which seeks to maximize the underutilized land for agriculture. Finally, scenario four represents an extreme vision in which the entire agricultural Frontier is occupied even if it is a protected area. Below, a brief description of each of the scenarios in terms of percentage change in the Orinoco basin area of Colombia: Scenario-1: This "business as usual" scenario defined expansion of 1.51% of oil palm, 0.12% of rice, 0.08% of forestry and 0.04% of soybeans.
Scenario-2: This scenario established 14.27% of the landscape to become agriculture dividing the area equally between the five basic products (oil palm, rice, forestry, soybeans and pastures for livestock).
Scenario-3: This scenario proposes expansion of 10.61% of agriculture (which were divided equally between rice and soybeans), 13.53% of livestock, 6.97% of forestry and 12.54% of agroforestry crops (assigned to oil palm).
Scenario-4: This scenario proposes 53.68% of expansion of all productive activity (oil palm, rice, forestry, soybeans and pastures for livestock).
These scenarios were simulated using the hydrological model developed for this study.

Indicators
One of the limiting factors of agricultural development is water (Abunnour et al., 2016;Salman et al., 2020). Long periods of  Trujillo, 2018). During times of drought, the pressures on the water resources intensify and become more evident.
In this research, we carried out a comparative analysis of the scenarios using 95% percentile of flow duration curve (Q95) as an indicator, which is associated with the flows present in times of low water. The flow duration curve (FDC) synthesizes information on the distribution and characteristics of the flow at a site, providing information for optimally managing the river's water resources. In hydrology, a FDC is a function that describes the variability of flow at a specific site during a period of interest. Q95 represents the flows equaled or exceeded 95 percent of the time, accounting for the base flows, rather than those due to extreme rainfall events (Ridolfi et al., 2020). One of the many applications of the FDC is to quantify the capacity to satisfy catchment requests (Dingman, 1981). The idea is to compare the percentage changes in Q95 between the baseline condition (historical condition) and the four proposed scenarios:  where ΔQ j is the normalized percentage variation of the flows of the flow duration curve, associated with the 95% percentile between the baseline scenario and k scenario, Q 95,Base Line is the discharge of the flow duration curve associated with the 95% percentile of the baseline scenario and Q 95, Scenario−k is the discharge of the flow duration curve associated with the 95% percentile of the k scenario.

Performance of the Hydrological Model
The model exhibited satisfactory results in the calibration process ( Table 2). In more than 85% of the stations selected as calibration points, the model presented NSE values higher than 0.6 ( Table 2), which is acceptable according to the criteria established by Moriasi et al. (2007).
As an example, Figure 7 shows the graphic results of the calibration at the control point established by El Barro station [35017040]. It shows the ability of the model to reproduce the dynamics of the hydrological regime draining to this control point. It also expresses the performance of the model in estimating the flow values in each of the seasons of the hydrological regime. This performance is reflected in the coefficient of determination (0.89) which highlights the capacity of the model in terms of the representation of the hydrological regime (1 being the highest value that this metric can reach). Metrics such as PDIFF (63.1 m 3 s −1 ), the R4MS4E FIGURE 7 | Comparative graphs in the calibration (1) and validation (2) at El Barro station [35017040].(A) Dotty plot graph of the variation of the objective function in the evaluation range of the parameters (a is runoff, b is upper limit of evapotranspiration and soil moisture, c. is groundwater fraction, d is reciprocal of groundwater residence time, Q thresholdj and V thresholdj indicate the percentage of contribution in each of the directions of the interaction, T rfj and T frj indicates the percentage of contribution in each of the directions of the river-floodplains and floodplains-river interaction and Ext is the percentage of water demand that is supplied by groundwater).(B) Graphs of parametric regional sensitivity, in red the objective functions with the lowest performance of the model, in blue the objective functions with the highest performance of the model.(C) Uncertainty using the Generalized Likelihood Uncertainty Estimation method (Beven & Binley, 1992 The quality and quantity of the information also represent an important factor in the performance of the model. Poor information leads to models that spread errors, generating greater uncertainty in the results. The scatter plot presented in Figure 8A shows that the model has a low identifiability of its parameters, except parameters a and b that show a slight concavity indicating that they can be identified. The regional sensitivity graphs ( Figure 8B) present the cumulative distributions of the evaluated parameters ranked from best to worst in terms of the objective function. With the GLUE methodology (Beven and Binley, 1992), the evaluated parameters are divided into ten groups of equal size according to the value of the objective function. The cumulative distribution of each group is represented as the probability value compared to the parameter values. The higher the gradient, the higher the sensitivity of the parameter. For this particular case, the only parameters that are sensitive are a and b. Figure 8C shows that the 90% confidence band contains most of the observed data, indicating that the structure of the model is capable of correctly representing the measured data. However, there are observations in periods of high flow, which are outside the confidence limits. This situation may be due to the uncertainty that exists in the measurement of flow rates, since the estimation is already made by means of a rating curve (Sikorska and Renard, 2017), a characteristic widely evidenced in hydrological modeling exercises (Ibarra et al., 2016;Shen et al., 2012;Yang et al., 2018). Another factor that can explain this behavior is the incomplete representation of the spatial variability of precipitation, as a consequence of limited monitoring in most of the Orinoco basin in Colombia. Figure 8C shows that the confidence band is wide, which denotes the parametric uncertainty that the model presents. Figure 9 shows the spatial distribution of the normalized percentage variations of Q95 for the four land use transformation scenarios.

Scenarios
These results show that the changes in coverage and land use derived from the four scenarios generate an important impact on Frontiers in Environmental Science | www.frontiersin.org October 2021 | Volume 9 | Article 673215 the water availability of the Orinoquia (with reductions greater than 50% in discharge during dry periods). In order of impact, scenario one turns out to be the mildest, while scenario four is the most drastic.
In scenario 1, most of the agro-industrial expansion is concentrated in a narrow northeast to southwest band of the flooded savannas of the departments of Casanare and Arauca, causing reductions of more than 50% in discharge during dry periods (shown in red in Figure 9, Scenario-1). However, in departments like Vichada, where protected areas and natural ecosystems predominate, the results showed increases in Q95 ranging from 5 to 25% in most HAUs. At the regional level, this scenario causes a 28% reduction in Q95 in the Meta River at the mouth with the Orinoco river.
The results of scenario two showed that the impacts on the water resource are greater than in scenario 1, extending throughout the department of Arauca, Casanare, part of Vichada and Meta. In this scenario, the Meta River at its mouth with the Orinoco River shows a reduction of 83% in Q95, in addition, tributaries such as the Guaviare River (at the confluence with the Inírida River) and the Vichada River (at the confluence with the Orinoco River) present reductions of 10 and 15% in Q95 respectively.
Scenarios three and four are the most unfavorable in terms of water resources given their expansion objective of 151,528 km 2 and 186,380 km 2 respectively. These two scenarios show a reduction of more than 50% in Q95, in 50% of the entire Orinoco basin in Colombia. In scenario 3, the Meta and Vichada rivers at their mouths with the Orinoco river and the Guaviare river at the confluence with the Inírida river, show reductions in flow of 95, 98 and 50% in Q95 respectively, while the same rivers for scenario four show reductions of 97, 95 and 48% in Q95 respectively.
An important characteristic that is general to all the evaluated scenarios is the extensive development evidenced over the flooded savannas, mainly by livestock and rice and oil palm crops. These developments affect the connectivity and water regulation capacity of the watersheds. They result from extensive logging, removal of vegetation cover, soil plowing and burning processes, all of which reduce the forests of the foothills and plains (Svenson, 1996;Tamaris et al., 2017), Although scenarios three and four present the greatest effects, the soil's physical and chemical limitations make these scenarios very costly economically, but not necessarily less likely (BuriticáMora et al., 2015;Osorio-Peláez, 2015;Bustamante, 2019).
The results obtained in the scenario analysis indicate changes in the hydrological regimes of the Colombian Orinoco basin that can generate direct impacts on biodiversity and ecosystems. The unique ecosystem of the flooded savannas support about 250 species of mammals, 567 species of fish, 761 species of birds, more than 700 species of plants, 108 species of amphibians and 119 species of reptiles in the Colombian Orinoco basin (Mora et al., 2015;Buriticá, 2016). A potential agroindustrial expansion in the Orinoquia region can reduce flooded savannas inundation and could have a great impact on the density and size of the populations and organisms that depend on flooded savannas (Kingsford and Thomas, 1995). For example, during 2014, in the department of Casanare, a region where more than 50% of its area is flooded savannas, there was a decrease in the quality of the habitat for Hydrochoerus hydrochaeris (capybara in English and Chigüiro in Spanish, the world's largest rodent) due to the scarcity of water, which generated a physiological pressure on the population and the subsequent death of a large number of individuals (Semana, 2014). The availability of water is one of the most important factors in the quality of habitat of the Hydrochoerus hydrochaeris (López-Arévalo et al., 2014). An agroindustrial expansion without considering the ecological limits of sustainability for ecosystems and species may generate great impacts that could put this delicate ecosystem at risk.
New water demands can generate reductions in the magnitude of the base flows, leading to biodiversity and ecosystem change. These water availability changes might include a decrease in the size of a population or local extinction of species (Kupferberg et al., 2012), reduction of species richness (Benejam et al., 2010) or terrestrial species invasion in the riparian and canal zone (Stromberg et al., 2007). Due to limits on water availability, agriculture and biodiversity conservation may be in conflict, especially during dry periods in areas with large monocultures. Poor water quality is another factor limiting usefulness of water. In rural areas of the Orinoquia, the population collects water directly from the rivers near their homes. Excessive use of agrochemicals in rice and palm crops may limit domestic water use (Reyes, 2020). This same situation affects ranchers who practice extensive cattle ranching.
The hydrological model demonstrated in this study allows decision makers to assess in advance potential conflicts related to the quantity and availability of water. It provides information on the amount of water in the 4,273 HAU, data that is essential for water allocation and planning. The model also allows us to understand how flows are reduced month by month when the water requirements of the hydro dependent water use sectors are increased, not only in the units where water is extracted, but also downstream of these units. If the model is combined with environmental or ecological analysis, planners can more efficiently allocate water across the basin and throughout the year, minimizing negative effects on ecosystems.

CONCLUSION
The coupling of the Thomas (1981) model and the Angarita et al. (2018) floodplain river interaction model demonstrated satisfactory results according to the performance metrics obtained in the calibration (NSE from 0.62 to 0.92) and validation (NSE from 0.49 to 0.92) exercises. Although the results of the sensitivity analysis showed that the parameters of the river-floodplain interaction model are not sensitive, the model allows, in a simplified way, to represent the natural dynamics that occur in the Colombian Orinoquia. In addition, the results of the Generalized Likelihood Uncertainty Estimation (GLUE) demonstrated that the model can structurally represent the flows recorded in the Colombian Orinoquia.
Sensitivity and parametric uncertainty analyses are necessary for any valid assessment of hydrological processes. Tools such as MCAT facilitate these processes, improving our understanding of how the model reacts to input signals and if its equations are adequate to represent hydrological records in the field. The adopted calibration scheme allowed the model to be parameterized in such a way that the representation of the hydrological dynamics of the Colombian Orinoquia was consistent with the records in the hydrological stations used in the calibration exercise. However, despite the parameters being clearly effective for calibration, they do not fully capture the characteristics of the basin. This is directly linked to parametric equifinality problems. For future work, a calibration scheme could be explored that relaxes the requirements of the basin hierarchy scheme use in this study, in favor of orienting the model towards biophysical processes.
Future efforts building on our hydrological model could be focused on capturing physical properties and their interactions with management practices. Moreover, among the opportunities to improve the hydrological model is the need to implicitly couple the demands of surface water with the equations that describe soil moisture decay, as well as groundwater demands with equations describing groundwater loss. The incorporation of processes such as rain interception could also improve the representation of the hydrological cycle. However, the objective of maintaining parametric parsimony and simplicity in the configuration of the model should not be overlooked.
The Colombian Orinoquia is considered the last agricultural Frontier of the country. Unplanned and intensive development, where the objectives of agricultural production and those of the ecosystems are unbalanced, can have serious repercussions on the water resource (decreases of more than 50% in discharge during dry periods). Lack of planning could compromise water security and put the ecological integrity of the Orinoquia at risk, as evidenced by model results of scenarios three and 4. Therefore, it is essential to integrate hydrological analysis in planning processes. In addition, planning must consider multiple scales, with analysis at the macro-region level (such as those addressed in this study) and at more detailed scales.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
Conceptualization, JN, CR, TW; Investigation, JN, CR; Methodology, JN, CR; Writing-original draft, JN; Writingreview and editing, JN, CR, TW. All authors have read and agreed to the published version of the manuscript.