Extended Project

Figure 2 Coello basin model Source: [6]

Validation method

The validation process was done through a comparison of the water volumes available in a period of time simulated by each model (WEAP & Python). This comparison was verified through the method of the Root Mean Square Error (RMS). Therefore, for the validation of the simulation, it was selected a reference scenario, this reference scenario does not involve restrictions of reservoir storage or climate change, and it works with the actual data on the Coello basin.

In the second part of the simulation there was selected one more scenario, worked by Meiline on her WEAP model, which will be describe bellow:

• Current situation (Scenario Reference 1980-2013): this scenario uses the actual data.

• Impact of climate change and socio-economics changes (Scenario All change 1981- 2013): this scenario is adapted from current situation scenario with several changes as listed in the following:

1. Increasing of potential crop evapotranspiration

2. Increasing water demand

3. Changing head flow by using the most extreme projection

4. Increasing domestic and industrial factor as the changing projection of socio

economic.

The input data used for the python model was proportionated in CVS files, which contains the information below:

• A daily series of flow in cubic meters per second for the reference scenario.

• A daily series of flow in cubic meter per second for the All change scenario.

• A monthly series of irrigation demand in cubic meters per second for references and all change scenario.

Case of Study

Coello River is located in the southwest of Bogota City, flows through Tolima Department and ends in the Magdalena River. The Coello Basin has several sub basin, they are Rio Toche, Rio Coello, Rio Combeina, Rios Barmelon, Rio Cocora and Quebrada Gallego. Water in Coello basin is mainly used for industrial, agricultural, livestock and public water demand [8]. Ibague as the capital of Tolima Department relies its water demand for domestic and industrial on the Coello basin.

For development of the model were used two CVS files provided by the WEAP model of Meiline Siahaan. One of those CVS files had a daily time series of flow in cubic meter per second and the second one had a monthly requirement of water for irrigation, both file belongs to the Coello basin in Colombia and were obtained from a meteorological IDEAM’s stations.

Figure 4. Water available at the exit of Coello basin. Scenario Reference a) WEAP model b) Python model

Figure 5. Water available at the exit of Coello basin. Scenario All change a) WEAP model b) Python model

Conclusions

The model developed in python managed to recreate the results obtained with WEAP models in a satisfactory way since the errors calculated are below 2%. This means that the program could obtained similar water allocation process than WEAP simulating the same behavior of the water demand, supply delivered, unmet demand and above all the behavior of the hydrograms at the end of each branch and at the end of the chosen basin.

This model does not have the GIS functionality, therefore the user must enter the data manually and check the mass balance equations for each branch. After that the user will be able to check the same results that was show before.

Additionally, within the process of water distribution in WEAP, it was observed that once this program makes the mass balance, if at some point it fails to supply the demand completely, it assumes that the precipitation was taken in a timely manner. Therefore, this program has continued its development, improving the spatial distribution of rainfall through current rainwater runoff transformation models, and also improving functionalities such as the simulation of the irrigation demand and the GIS functionalities.

Figure 2. Cells of the Global Climate Model: MPI-ESM-MR, in the Bogotá river basin. Coordinate reference system: WGS-84

For the local data, 47 ground rainfall gauges stations selected from IDEAM (Institute of Hydrology, Meteorology and Environmental Studies, Colombia) were used for the period 1958-2016, as shown in the Figure 3, observing that there is an average of 12239 daily precipitation data per station, for use as predictors, with an average value of 2.76 mm, an average standard deviation of 6.69 mm and an average maximum value of 105.54 mm.

Figure 3. Ground rainfall gauge stations selected for the downscaling in the Bogotá River Basin

The Chaotic Statistical Downscaling model (CSD) evaluate the presence of deterministic chaos for different rainfall accumulation intervals for both of the time series (rainfall gauge stations and the global climate model). Then a chaotic predictive model is constructed with the results of the time delay, the embedding dimension and the Lyapunov exponents found for the "optimal" precipitation accumulation interval of each dynamic system. The predictive model is based on the synchronization between two dynamical systems, given by the parameter μ of the mutual false nearest neighbor’s method (MFNN).

According to (Sivakumar & Ronny, 2010): "The phase space is essentially a graph, whose coordinates represent the variables necessary to fully describe the state of the system at any time, the trajectory of the phase space diagram describes the evolution of the system for an initial state. The "region of attraction" of these trajectories in the phase space provides important qualitative information when determining the degree of complexity of the system ", as shown in Figure 4 where the phase space is observed in three dimensions for different accumulation intervals of precipitation. The limiting set that brings together asymptotic trajectories close to equilibrium is known as 'attractor'. The attractors of deterministic chaotic systems can exhibit an unusual type of self-similarity and present structures at all their scales and it is therefore necessary to find an appropriate dimension of the phase plane, such that the structure of the attractor remains invariant.

Figure 4. Comparison of the phase space of the precipitation time series for Cell 1 of the Bogotá river basin for rainfall accumulations of 1 day (left), 5 days (center) and 30 days (right).

For the evaluation of the presence of deterministic chaos in the dynamical systems Lyapunov exponents were used, for which it was necessary to reconstruct the phase space with the Method of Time-Delay, which finds the appropriate values of the time delay (τ) and embedding dimension (m) to capture the attractor dynamics. The autocorrelation function and the mutual information were used for the selection of the time delay, and the embedding dimension was selected using the correlation dimension method, the False Nearest Neighbors method (FNN) and Cao’s method, which were successfully used in (Gallego, 2010), (Hernández, 2009) and (Siek, 2011).

Within the different phenomena of dynamical non-linear systems, the synchronization of dynamical systems has had a great importance in the last decade, being used in a wide variety of practical and experimental applications, such as electrical and electronic circuits, lasers, telecommunications systems and chemical reactions. The synchronization of two or more dynamical systems is a fundamental phenomenon that is obtained when at least one of the systems changes its trajectory due to a coupling (unidirectional or bidirectional) with another system, allowing a coherent behavior in coupled systems. In chaotic systems, given the dependence of the initial conditions and their evolution in different attractors and dimensions, the synchronization process allows the systems to follow a common and defined trajectory.

The mutual false nearest neighbor’s method created by (Rulkov, Sushchik, Tsmiring, & Abarbanel, 1995), is a statistical technique based on the calculation of the parameter μ, which evaluates the local neighborhoods between two time series, so that μ takes values of the order of 1 if exist complete general synchronization, otherwise, μ has to be a number whose magnitude is comparable with the product of the attractor size divided by the product of the distance between the nearest neighbors in the time series x and y. The formula for the calculation of μ is given in Equation 1 and the relation between points of the attractors is shown in Figure 5,

Equation 1. Mutual false Nearest Neighbors

Where,

Xn and Yn , are points of the drive and response system at an instant n

XnNND and YnNNR are the closest neighbors of and in their respective systems.

XnNNR and YnNND are points of the opposite system at time n for ynnR and xnnD.

Figure 5. Relation between Xn ,Yn ,XnNND , YnNNR, XnNNR and YnNND for drive and response systems

Results and Discussion

With the purpose of obtain a wider range of analysis, the evaluation of the presence of deterministic chaos was performed in intervals of 1,3,5, 7, 10,15 and 30 days. Analyzing the results, it is observed that in the cells of the GCM and in the rainfall gauge stations time series, starting from a rainfall accumulation of 5 days the type of movement of the dynamic system is no longer random (noise) and becomes mainly deterministic chaos, ensuring the short-term predictability for climate projections. It was also observed that in both time series the type of movement of the dynamic system changes from deterministic chaotic to stable for a rainfall accumulation interval of 15 to 30 days, ensuring the existence of long-term predictability for monthly climatic projections.

The phase space reconstruction for local gauge stations and GCM cells shows that the embedding dimensions are 6 and 8, respectively. This indicates that the GCM is more complex than the dynamic system shown by the local gauge stations. The complexity is also reflected when calculating the average attractor size of 3.5 for the local stations and 4.86 for the GCM. These results are of great importance, since they allow to identify the degree of difficulty to synchronize the two chaotic systems.

The values obtained for the parameter μ of the mutual false nearest neighbor’s method are greater than the theoretical value of μ = 1 ideal for general synchronization. However, it is possible to observe that using the CSD model it was possible to decrease the difference with the ideal theoretical value by finding the minimum possible distances for - and - , which represents a slight increase in the synchronization of the two systems, as presented in Figure 6 for the station pr_21185040

Figure 6. Results of the parameter μ of mutual false nearest neighbor’s method for station pr_21185040

In the process of compare the downscaling model, a wide variety of statistical downscaling methods were used in the Bogotá river basin, for a total of 16 techniques, including: CSD, k-NN Bootstrapping, Delta methods, Analog methods, Quantile Mapping method, Weather Type methods and Generalized Linear Methods, most of the methods were include in the Meteolab toolbox (Cofiño, et al., 2013). These methods were evaluated under three different measures of error (Mean value difference, Max value difference and Root-Mean Square Error (RMSE)), in the Figure 7, the RMSE comparison is shown.

Finally, it was possible to compare the monthly multiannual precipitation, the annual precipitation, and the return period for a wide variety of downcaling techniques with the reference period, and use these results to estimate the future river flows using a hydrological model. In the Figure 8 is shown the spatial variability for the percentage difference between the historical reference and the projected time series using the CSD model.

Figure 7. Comparison of the RMSE error obtained for different downscaling techniques

Figure 1. Conceptualization for TACD tank modeling.

Each cell in the model domain is classified in one of the previous categories, according to some criteria defined by the researcher. In the Brugga basin, criteria for classification is shown on Figure 1. As mentioned before, for each dominant process, capacity, conductivity and permeability parameters are calibrated to appropriately represent hydrological behavior of sub basins where water flow and water level is available for a convenient period of time.

Input of climatic data depends on the represented variable. While precipitation and temperature are spatial and temporally distributed data, land elevation and average time of solar radiation per day are spatially distributed data that remain constant through the model evaluation. Other data such as observed water flow and water level are scalar time series, that is, they are considered a boundary condition for a single cell on the watershed.

Study Area: Coello River Basin

Coello river basin is located at the center of the Tolima Department, in Colombia, with an extent of around 1800 km2, representing 7,8% of the total area of Tolima. Coello river starts at the Don Simón moor (3850 m.a.s.l.), flowing with an average slope of 3.2% through approx. 121 km from west to east until it reaches the left bank of the Magdalena River (300 m.a.s.l.).

Coello river and its tributaries play an important role in the regional development, because several nearby cities get water from it, such as Cajamarca, Ibagué, El Espinal and Flandes, among others. Its main tributary rivers are Toche, Bermellón, Cocora, Andes, Gallego and Combeima. Coello river has an estimated water supply for the region of 30,9 m3/s and an estimated water demand of 23,8 m3/s. Main water demand sources are drinking water supply, crop cultivation (coffee, grains, citric fruits, sugarcane, among others), animal husbandry (meat, eggs, eggs, etc), and mineral extraction.

Coello river basin has a mean annual precipitation of 1517,8 mm, on average being April the rainiest month (188,9 mm) and January the driest (66,7 mm). Mean temperature for the basin is 19.8°C, but varies from 28,6°C in the lowest points to 3°C at the moor. In modern times, flood events have been observed at the study region, such as those that occurred in Perales stream (1989), Anaime River (1993) and El Oso stream (1998), causing several material losses.

Figure 2. Coello basin location in Colombia, with main rivers and elevation model.

Methodology

Methodological Description

In order to use the TACD model to represent hydrological behavior of the Coello river basin, the model was updated to a recent version of its GIS software PCRaster version 4.1. Some capabilities of previous versions were no longer supported, and some new features were added. Therefore, a complete code rewriting in Python was necessary. In the current version of PCRaster, a Python framework interface is available for use.

Having in mind that rewriting the model scripts could lead to numerical differences between the two versions of the TACD model, validation of the new results with the original model results were necessary to ensure correctness of updated model logic, and for that reason, the Nash-Sutcliffe Efficiency was used between the original and the updated model in the Brugga basin.

Being sure that the updated model accurately calculates hydrological behavior for a basin, then the model and the input data for Coello river basin should be adapted to one another, such as maps and time series formats, folder organization, initial parameters, etc.

After setting up the Coello river basin model with TACD2, a calibration process followed in order to obtain a subset of model parameters that minimizes the error between observed and modeled runoff at the Payandé flow station. Results and comments on this process are described below.

TACD Model Main Features

As mentioned before, TACD is a distributed hydrological model that runs several routines for the hydrological processes in a watershed on a daily basis and outputs spatial as well as scalar information regarding soil moisture, snow content, stream flow and evapotranspiration, among others (Roser 2001). It is important to remember that as a distributed model, TACD calculates each of these routines averaging for each cell (pixel) in the study area by using several kinds of input such as spatio-temporal data (precipitation, temperature, vegetal coverage maps, Leaf-Area Index), spatial stationary data (channel width, digital elevation model, forest areas, water flow direction, manning’s n coefficient, streams, sun hours per day, urban areas) and scalar data (observed waterflow, tank parameters), as required in each routine:

• For the potential evapotranspiration, the Turc-Wendling approximation is used, which takes as input the daily solar radiation and the average daily temperature for each cell in the study area

• For the snow content balance, a 2-tank submodel was developed. This submodel separates the solid and the liquid part of the precipitation and calculates freezing and melting processes

• For the interception of precipitation (mostly due to plant canopies), a tank submodel is used. In this submodel, a fraction of the flux going towards the ground is intercepted in a small tank with limited capacity, so that if the water volume plus the current content of the tank exceed the tank capacity, then surplus will continue its way to the ground, as if not intercepted. Otherwise, evaporation processes occur on the tank and eventually might empty the tank. For this submodel, Leaf-Area Index (LAI) and Normalized Difference Vegetation Indices (COV) maps are required

• For the runoff in urban zones, the precipitation has a short concentration time compared to natural fields and therefore the model concentrates directly that precipitation into the most upstream part of the closest river where this urban subbasin drains to

• In saturated zones (e.g. lakes), infiltration rates are low compared to runoff and therefore, most precipitation actually flows directly with little retention effect. This is achieved with a two tank submodel whose upper tank capacity is low and whose lower tank still allows subsuperficial flow. In a similar way, the part of the precipitation that falls on the stream surface, is directly added to the waterflow in the stream

• For the soil infiltration, routines are very similar to the HBV model and represent infiltration as an exponential function of soil moisture and field capacity. In this routine, actual evapotranspiration is also calculated as a function of soil moisture, the field capacity and the potential evapotranspiration

• For the runoff routine, the tank parameters are used in order to calculate for each cell how the water flows through each of the tanks in each cell, according to the runoff-rainfall transformation defined for each cell

• For the channel flow routine, water routing is calculated along the river. This is done through the kinematic wave implementation in PCRaster Python Framework, taking as input the manning’s n coefficients, channels width and average length of channel per cell, as well as the spatial and temporal discretization dx and dt

TACD2: Model Update to Python

When writing the Python version of the model, some features of the original PCRaster version were adapted with built-in functions of the Python language. Besides, as the new framework is object-oriented, some major changes in model structure were made for readability and completeness.

A class called TACD2 is instatiated by setting up the path to the clonemap (study area) in .map format. After that, a DynamicFramework must be created for the TACD2 object with the initial and final timestep (with Python arguments firstTimeStep and lastTimeStep), and then the dynamic model is run by calling its run() procedure. All other procedures and attributes inherited from the PCRaster DynamicModel class are also available, such as setQuiet() for enabling or disabling output for debugging.

Figure 3. TACD2 script structure (object-oriented).

Model Validation

For testing goodness of fit from the updated model to the previous model, Nash-Sutcliffe model efficiency coefficient (NSE) was used as described in Nash et al. (1970). NSE, calculated as in Eq. [1], can take values from –∞ to +1, where an efficiency of 1 (NSE = +1) corresponds to a perfect match between observed and modeled data. Moreover, an efficiency of 0 corresponds to data where the modeled data is as precise as the mean of the observed series in forecasting the observed results.

Where:

𝑄𝑡𝑚: Modeled variable through time t (e.g. flow)

𝑄𝑡𝑜: Observed variable through time t (e.g. flow)

𝑇: Total time, i.e., number of timesteps

𝑄𝑜̅̅̅̅: Mean for observed variable

Model Setup for Coello River Basin

Taking into consideration difficulties by updating TACD2 model, a new folder structure for the model was defined, as shown in Figure 4. Once correctness of the TACD2 model was validated, spatial and temporal information for the Coello basin was to be converted to PCRaster required formats.

Other spatial information required for Coello basin was urban settlements location, soil classification, Leaf-Area Indices (LAI) and Normalized Difference Vegetation Indices (COV). Spatio-temporal data such as precipitation and average temperature was gathered from the European Center for Medium-Range Weather forecasts (ECMWF), ERA-Interim Database in NetCDF format. Time series such as water flow and water level in water level stations of study area were gathered from the Colombian Institute of Hydrology, Meteorology and Environmental Studies IDEAM, through its information service Surface Waters Observatory.

Figure 4. Folder structure of TACD2 model.

Parameter Calibration

Calibration scheme of TACD2 to Coello basin took into consideration water flow and water level at the Payandé Station. During calibration process of parameters, CPU time needed for this task was a bottleneck. A total of 67 parameters have to be adjusted. Most of these parameters are related to runoff generation, that is, tank model constants (Tank capacity, lateral flow between tanks, infiltration and percolation). As an example of the calculation complexity, the chosen time series take approximately 1 hour for each model run.

In order to optimize the calibration process, the standard python library scipy was used, more specifically the standard minimize procedure was applied to the Root-Mean Squared Error (RMSE) between modeled and observed waterflow at the Payandé Station.

Figure 5. Comparison of intercepted precipitation series.

Figure 6. Comparison of intercepted evapotranspiration series.

Figure 7. Comparison of cumulative intercepted evapotranspiration series.

Model Setup and Calibration

The first step for validating the calibration, the comparison of observed and modeled flow at the Payandé station is verified. For these series, the NSE coefficient is -0.57296 (RMSE=22.4937). However, Figure 8 shows that the mean trend of the modeled results corresponds to the mean trend of the observed data, but variability in both series is very different. A longer calibration process is encouraged in order to achieve a better fit.

Figure 9 shows intercepted evapotranspiration (blue) as well as intercepted precipitation (orange). The intercepted precipitation is the sum of the precipitation that enters into the interception model tank and the interception surplus of the previous evaporation step. From the figure, it is observed that every drop of water that is intercepted evaporates and goes back to the atmosphere. Further calibration of the parameters involved in this process is necessary.

Figure 9. Intercepted evapotranspiration and intercepted precipitation.

Figure 10 shows water content in soil for the final timestep of the TACD2 model (left) and the soil classification for the basin (right). Problems due to spatial resolution of temperature can be seen in this figure, because the grid resolution for this map dataset was 11km x 11 km, compared to the 100mx100m resolution of the elevation model. Besides, a big contrast in the upper part of the basin can be observed. This water content difference is caused mainly by the temperature difference between adjacent cells. Hydrological classifications also play an important role in the model, as it can be observed that some zones become “highlighted” in the water content due to change in the soil classification.

Figure 11 displays final results for the modeled water flow at each station in the study area. Although only Payandé Station was used for calibration, it can be seen that most discharges from other subbasins have similar shapes at different scales with some small variability to each other as an expected result of local processes at each subbasin. From the current modeled flow at each station, it can be shown that the hydrological water supply potential could be estimated for each subbasin as follows: At the Payandé Station, 67.4 m3/s; at El Carmen, 42.8 m3/s; at Puente Carretera, 33.0 m3/s; at Yuldaima, 15,7 m3/s; at Bocatoma (intake), 0.7 m3/s; at Puente la Bolivar, 3.8 m3/s; at Puente Luisa, 2.1 m3/s.

Figure 10. Soil water content for last timestep (up) and hydrological classifications for Coello basin (down).

Figure 11. Modeled water flow at each station information Coello river basin.

Figure 12 shows water balance at the Payandé Station through time. One can see that there is an almost steady decrease in water balance as time advances, up to a water balance difference of -2121,38 mm at the end of the year. When comparing it to other stations, such as Puente La Bolívar or Puente Carretera, all of the water balances show the same decreasing trend, probably caused by the initial values of storage, which were calibrated for the starting part of the simulation runtime but may also be calibrated to ensure a zero sum in water balance.

Figure 12. Modeled water flow at each station information Coello river basin.