1. Introduction
The Bogotá River Basin has a total area of 276.2 km2, which is constituted by a unimodal rainfall regimen between April and October. Along its route, ecosystems such as Andean Woods resulted being vulnerable due to water contamination. There is a compelling need to recover and to conserve this resource [1]. Tanneries are industries belonging to a traditional economic activity in Cundinamarca - Colombia, where animal raw hides are transformed in leather. According to the POT (Plan de Ordenamiento Territorial), up to 131 tanneries have been established around the riverside from 1,820 to 2,000. Most of them came from Villapinzón, Cundinamarca. Nevertheless, the number of small tanning companies started to decrease from 2001 because of closes and sanctions performed by the Corporación Autónoma Regional (CAR), which monitors and regulates environmental policies of their discharges [2].
Water harvesting and discharges made by tanneries in the High Bogotá River Basin have been invariant due to the significant volume needed by these industries, which is approximately to 12.30 m3/t [2]. According to Observatorio Regional Ambiental y de Desarrollo Sostenible del Río Bogotá (Orarbo), the High Bogotá River Basin has a good water quality index (0.90) in upstream of Villapinzón and regular index (0.37) in downstream of Chocontá [3]. However, the main issues are that Cr III composes these discharges and there is an absence of adequate treatment by 98% of the industry, which are conformed mainly by small companies [4, 5]. The risk that hosts a high concentration of Cr III in the water is because it can react to Cr VI, which is highly harmful and is classified as a carcinogenic substance [4]. This situation in public health is important and trends to be economically important; thus, regulations and management should be implemented [5].
Simulation of persistent pollutants in Colombia has not a broad background. Preview studies are restricted to water bodies in high mountains, and only few have been developed. The software Water Quality Analysis Simulation Program (WASP) 7.0, promoted by US Environmental Protection Agency, predicts and performs simulations of water quality in response to natural phenomena and human-induced pollution. Nevertheless, this program has not yet been applied in Colombia. Simulations in the Suquiapa River in El Salvador and Tibacay River in Ecuador demonstrate an accurate adaptation of WASP to segments, hydraulic variables, and boundary conditions from water bodies in high mountains, such as High Bogotá River Basin [6, 7] WASP includes submodules which allow studying different compounds, such as the simulation of Cr III through simple toxicant submodule [8].
This study focuses on the simulation of Cr III, which is a persistent pollutant and a non-biodegradable mineral harmful for human and animal health. Hydrologic dynamics of Cr III provides monitoring and control strategies for sanity and environmental impact reduction of heavy metals. This simulation in the High Bogotá River Basin affords to understand changes in Cr III according to spatial and temporary factors. Therefore, the aim of this research was to elucidate these changes along this water body via WASP simulation. An initial validation was developed by both pollutant load error and continuity of water volume error. A quantitative validation of the Cr III concentrations was performed using the Index of Agreement (IOA) validation parameter.
2. Materials and methods
2.1 Delimitation and characterization of the study zone
This research was developed for the High Bogotá River Basin, which is located in the north-eastern sector of the Department of Cundinamarca [Figure 1]. 13 segments, 7 tanneries and 2 streams of discharge were assigned [Figure 2]. Values of water level were extracted from level control and CAR monitoring stations from El Triunfo and Villapinzón.
The High Bogotá River Basin is constituted by the municipalities of Suesca, Lenguazaque, Machetá, Villapinzóm, and Chocontá. In 1,982, the source of the Bogotá River was declared a forest reserve zone, reaching 1,128 ha. Its headwaters are located about 3,450 m asl in the wasteland Guacheneque in Villapinzón, and its mouth is about 2,600 m asl in the Sisga River. The median annual precipitation is 850 mm, being July and January the rainiest and driest months respectively. The average annual relative humidity has a similar behavior of the precipitation, following a monomodal regime. The highest relative humidity values were registered between April and September, with July being the wettest month with a value of 91% [1].
2.2 Selection of Water Quality Analysis Simulation Program (WASP)
In this study, WASP ® was chosen to model the contamination of Cr III in the High Bogotá River Basin. This selection was based on the comparison of license type (commercial or free), model type, water body admitted, and the phenomenological variables included of available programs. According to previous studies in rivers with similar characteristics [8, 9], the user-friendliness on the interface, and the simulation scope for conservative pollutants were selected as essential factors in this study case. Also, its easiness to install in any computer made the selection easier and remarkable for future investigations.
2.3 Model and equations
A partial differential equation based on the law of mass conservation was formulated including terms that describe the pollutant transport, inputs and outputs of the system, variables related with hydraulics conditions, diffuse discharges, punctual discharges, and the analyte capacity (Equation (1)]. This equation presents the continuity of the fluid in a Cartesian coordinate 3D system and make reference to the lateral (x), longitudinal (y) and vertical (z) transport, which can be produced by the mechanisms of advection (uc, vc, and wd), diffusion (E), loads direct or diffuse (S L ), exchange of matter with sediments (S B ), or kinetic transformations (S k ).
An equation was developed from variables associated with X-axis and time (t) (Equation (2), Table 1]. This equation resulted from the suppression of the terms related to Y- and Z-axis, which comes from a vertical and horizontal homogeneity assumption in each segment [9] WASP® applies a numerical approximation using finite differences, with centered numeric derivatives, in first-order Eulerian scheme, for the solving of the mentioned equations [9]. This numerical approach is based on an Euler method for Ordinary Differential Equations (ODE). This approach is numerically implicit and is the basis for the finite differences scheme used by the model.
The pollutant transport in a water body is dependent on the dispersion coefficient (Et). According to Barchiesi and co-authors (2015), the Fischer equation (Equation (3)] provides the best estimation for the current study case [11].
Flows were determined from the Manning equation (Equation (4) ). Variables such as the geometry of the section [Figure 3], the wet area (Equation (5)], the wet perimeter (Equation (6)], and the hydraulic ratio (Equation (7)] were taken from Villapinzón-2120815 limnometric station and the classification proposed by other authors [Table 2] [12].
The parabolic sections scheme shown in Figure 3 was used to model the system based on Equations (5), (6), and (7). This geometry section is based on the Schumm classification [13].
The Peclet number was calculated (Equation (8)] [14], and the dominating mechanism was determined in the water body taking account guidelines from other studies [15, 16] criteria. These criteria established that an advective runoff trends to present an infinite Peclet’ number and advective dominates over diffusion starts when Pe > 1. Therefore, a highly convective flow was identified in the studied system.
2.4 Discharges concentration
The concentration of the discharges was determined from the stochiometric adjustment of the chemical reaction of Cr in the tanning process.
From the balanced reaction (Equation (9)], the concentration of Cr III in the discharges from the tanneries was estimated from the mass percentage of Cr in the products, the water consumption reported in the balance performed by US EPA in 2006 [15], and the number of processed hides [16]. The output flowrate of the tanning process is 2.36111x10-6 m3/s of process water, 0.221 kg of chromium sulfate are commonly used, with this data and the mass percentage of Cr III (3.32%) present in the used salt during the production process, the concentration of Cr III in the stream is obtained (6890000 µg Cr/L) [17].
2.5 Validation
The error percentage in the simulation was determined by the calculation of continuity error for the water volume and pollutant in the segments in the simulated stretch. The continuity error in the volume corresponds to the difference in the initial water mass in the simulated stretch (S-1) with the sum of the flows in the inputs from the system (D-1 to D-7), including the value of flow in the final segment (S-13) and San Pedro’ and Chingacio’ streams (T-16 and T-20, respectively) (Equation (10)]. The error of the pollutant charge in the mass of the system was calculated based on the concentration of Cr in the system’ discharges which was an input to the model (Equation (11)]. Chrome prediction results at the end of each segment were subtracted from the concentration of the discharges, and then, divided by obtained value in the simulation. The Index of Agreement (IOA>0.6) was applied to evaluate the performance of the model [18]. Equations (10) to (13) illustrate the error calculation and validation parameters considered for the model performance validation.
The data used for the validation consists in the observed and simulated Cr III concentration (Equation (11)], data used is available through the link in the supplementary material of this article.
2.6 Application of WASP and configuration of the model
The simulation was developed using WASP model v 7.0 [19]. The simulation was performed on a computer with 6 GB RAM, Windows 10 as the operating system, and Intel Core i5 - 5500U. The model set-up was based on the Net flow option, this option assumes one single flow direction for the water in the simulated systems [19]. Euler method was selected as numerical solution method, it allows solving ODE along the segments distance considering specific initial and boundary conditions [20]. Additionally, the Euler method presents a fast and good convergence. The simulation execution time was 24 hours using this numerical approach.
The sub-model simple toxicant was chosen for the system physics representation since it has affinity for metal compounds modeling, among other pollutants [6]. Depth, width, length, Cr concentration, and segment type were provided to WASP as input information for the simulation. The volume and characterization of each segment were based on the measured data collected data.
Advective and dispersive factors were defined as follow: advective transport were established by entering the flow values at the initial state of each segment allowing to represent the flows associated with the water streams inside the simulated system. Dispersive factor was determined based on the cross and longitudinal dimensions. The post-processing module of WASP allowed to export and analyze concentration changes of the Cr III at spatial and temporal scales [Figure 4].
3. Results and discussion
The concentration of Cr III and flow values associated with discharges points [Table 3] were obtained from the empirical equations and experimental data [15]. The concentration used in Equation (11) was calculated by multiplying flows and concentrations of discharges. The error percentage for pollutant load was minimum with the model implemented, and accurate results were obtained based on the relation between the simulated and observed concentration. An invariant pollutant load was calculated based on the assumption that all discharges have a Cr-III-rich tanning origin. Pollutant load values are agreed with observed pollutant concentration (6890000 µg Cr/L, 2,36111 E-06 m3/s) [15].
Table 3
Concentrations and time where a maximum dilution (mixing condition) of Cr concentration was achieved
Figure 5
Cr concentration, flowrate profiles and error percentages of Cr along the simulated stretch of Bogota River
Figure 5 shows the observed and simulated Cr concentration profiles, the water flowrate profile and the punctual percentage of error is reported. The errors computed are in a range between 0.14% and 39.26%. A correspondence in terms of trend for the observed and simulated Cr concentrations is evidenced, even so, the simulated Cr concentrations are more sensitive to changes in the water flowrates. By example in S-9 a decrease in the flowrates of water is registered, but by comparing the Cr concentrations the simulated profile shows a peak of higher Value meanwhile the observed Cr concentration does not change strongly.
From Figure 5, it is possible to observe in sections S-6, S-7, S-8, S-10, S-11 and S-12 a good correspondence between the simulated and observed Cr concentration profiles. These sections of the river correspond to a small change in water flowrates.
Mixing zones were identified in each segment from the stretch. These zones represent the maximum possible dilution of Cr III in the water flowrate. Mixing zones also avoid fluctuations in the pollutant concentration. Mixing zones are constituted mainly by 2 steps: turbulence zone and far field zone. Turbulence zone is where the discharge is added, and the Far field zone is where changes in the pollutant only depend on the hydraulic phenomena. From the simulation, initial, minimal, and maximum concentrations of Cr were obtained as results [Figure 5 and 6).
A progressive increase of Cr concentration was observed in the relation of time until achieving the mixing condition in the system. From this behavior, an interruption was detected in S-6 were initial concentration, and the mixing zone reached values of 44,38 µg/L and 27,74 µg/L, respectively. This result is more likely due to the inlet of a stream (D-16) of 0.499 m3/s clean of Cr. This input dilutes Cr because the available water volume is greater. This situation is also evidenced in S-6, S-7 and S-8 for mixing condition.
S-9 presents a higher concentration at mixing conditions, near the initial concentration of the segment. This is presented because S-9 receives 2 discharges with Cr leading to an increase in the pollutant concentration. For S-10 and S-13, WASP did not achieve the stability condition, for this reason fluctuations of Cr concentration are represented by the model. We inferred that maximum dilution was not reached in the 24 h of running time for these sections (S-10 and 13).
Modification in the geometry of the river sections and flows of the discharge points were performed, variations in the Cr concentration were obtained for each segment. The advective mechanism is linked with the dispersion coefficient set in WASP with a value of 1 m2/s. A modeling study developed in the Tumaque basin in Venezuela [14] supports this value based on the domination of the advective transport in the system. The coefficients obtained with the Fischer equation (Equation (3).) have dispersed values for the system analyzed in all the 13 segments. Due to this, the Peclet number (Equation (8)] was computed, which makes it possible to identify accurately that the advective mechanism predominates in the transport of the contaminant [21, 22]. Detailed data is available through the link in the supplementary material of this article.
Variations in the flow and the initial geometry of the segments allowed establishing the sensibility of the Cr concentration to the model configuration parameters. This observation can serve as qualitative validation of the model performance. From Figure 5, an increase in the flowrates, mainly associated with river section depth, corresponds to a decrease in the Cr concentration along the stretch of the river. These results are coherent with the equation of continuity (Equation (2)].
Based on the obtained results [Figure 6], the hydraulic retention time (HRT) and the Cr concentration tend to increase along the stretch. 2 sections with local decrease were observed (S-6 and S-11), where pour Cr streams are added to the system. Profiles in Figure 6 allow observing the influence of the discharge flows and the geometry of the sections in HRT and concentrations. A clear relationship between concentration of Cr and HRT is also observed.
Cr concentration, along with segments from S-6 to S-8, was affected by the San Pedro stream. In contrast, concentration resulted higher along with segments from S-9 to S-11, which is explained by longer segment distances (837,994 m to 1083,835 m). An increase in the HRT affects Cr concentration because superficial absorption and sediment formation are favored, this behavior was reported in a previous work in the Bogotá river at Villapinzón sector where Cr traces were predominant in the sediment of the water stream [23].
Similarly, a validation process to establish the accuracy and reliability of the simulation performed was implemented based in 2 quantitative indicators computed between observed and simulated results of Cr concentrations and water flowrates: the error and the IOA. The error in the whole stretch was computed as 2.6 % average for the water flowrates [Table 4]. For pollutant concentrations, the error percentage in the validation was 13,19% average. These error values are low indicating a good representation of the measured data by the proposed model.
The 50% of the granulometry that composed the bed of the simulated stretch from the river is >2.7 mm [24], which indicates that infiltration can be avoided because evaporation or discharges are the main processes associated with loses in this system.
The error associated with the pollutant charge was higher in the segments S-6, S-7, and S-9 reaching values of 39,26%, 26,49%, and 19,29%, respectively [Figure 5]. The difference between concentrations and flows in the segments are the main cause for the increase of the accuracy percentage in the second validation. Previous modeling studies made with QUAL2K [10] showed NASH coefficient in its calibration that ranges from 0.02 (insufficient) to 0.9 (excellent) [25], which confirms some sections present high accuracy and others some variation in the water stream.
For the evaluation of WASP model performance, the IOA was computed as main validation parameter. The IOA value of 0.853 shown in was obtained from the simulated and observed concentrations of plotted in Figure 5. This IOA value indicates a high correspondence between observed and simulated data of Cr concentration in the water stream. This means WASP model implemented satisfies the criterion established by Willmont, reaching a greater value of 0.6 for this validation parameter [18]. IOA values were computed using OpenAir® Library in R programming language Figure 7.
4. Conclusions
WASP model was implemented in this research to simulate the Cr concentration dispersion at the High Bogotá River Basin. Results with 13% average of error for the pollutant in the water stream were obtained. The segments S-6, S-7, and S-9 presented the highest error regarding the pollutant concentrations with error values of 39.26%, 26.49%, and 10.29%, respectively. These error values are acceptable but can be improved for future studies by using highest spatial resolution and more simulated time. The model is capable to predict the hydraulic phenomena in the simulated river section with coherent results according to the dominating equations of the system. WASP model may get more accurate results by using hydraulic data from more accurate experimental methods and longer measuring campaign of the variables involved in the system. These experimental methods must be oriented to determine the geometry and water flow rates in the analysis system. The experimental method would get data capable to improve the WASP model prediction in terms of pollutant concentration error. Potential WASP model accuracy, also considering WASP model flexibility, make WASP to be a useful tool for the management of hydric resources based on mathematical and numerical modeling oriented to scenarios assessment.
The 13 segments used in this study present a gradual increase of the pollutant concentration of Cr. These increments are observed upstream the mixing zone. The last 4 segments did not fully reach the mixing condition presented a Cr concentration of 70.48 µg/L. This value corresponds to a high concentration since, in 2003, the Cr presence in surface water was between 1and 10 µg/L. This change puts in evidence a large impact by the industrial activity on the analyzed water system in this study during recent years.
The Cr III under high concentrations of organic matter can be oxidized to Cr VI [26], due to this feature Cr III can be taken as an indirect indicator of Cr VI a very harmful substance for human and animal health since it is a compound listed as a cancer promoter by WHO (World Health Organization). Environmental Control institutions in the region must develop sanitary regulations in the water system since currently, an optimum scenario for Cr oxidation is taking place in the river section. This situation is a potential risk and since Cr oxidation is mainly catalyzed by the organic matter heading to the river from sources such as residual domestic water and ledger industries residual water.
This study allowed validating the implementation of WASP model as a useful tool for the water resources modeling in Colombia and Latin America, due to the similarity in the hydric systems in the different Latin American countries, these systems are characterized by high mountains rivers or even lentic systems. WASP model versatility to perform different types of simulations with many variables allows the development of further important studies for water pollution. Considering the high rates of industrial increasing activities and water consumption in the region makes WASP has the potential to be important for future research, resources planning and development.
In this case, regarding the ledger production industry the model has significant opportunities to be implemented. The potential applications of the model cover use in the environmental assessment of polluted loads into the system and optimization of treatment processes for the industrial wastewater [26].
This research implemented a quantitative validation concluding WASP is capable of predicting accurately the behavior a system of water with polluted loads with an IOA of 0.853 between observed and simulated pollutant concentration [18]. The variables considered by the model as inputs and outputs can be used for assessment of impact in public health, sanity, and resource quality.
The modeling performed corresponds to a baseline which identifies errors and options of improvement oriented to a more accurate simulation in the future. This study is a starting point for related research in the topics of water resource quality assessment, scenarios analysis, and infrastructure development for decision-making processes in the region involving the effects of industry and local policies.
Supplementary Material
Table 5
Geolocation of segments, discharges and streams of the stretch of the river coordinates in UTM
Table 6
Available models for simulation of water systems quality
| Software | Advantages | Disadvantages | References |
|---|---|---|---|
| QUAL 2K | Useful for simulation of algal concentration and Biological Oxygen Demand in water. The discretization of segments can be handled irregularly. | Specific application for users with defined knowledge and specific skills. The model has limitations in simulation parameters. Heavy metals cannot be simulated. | [27] |
| QUASAR | Handling of persistent pollutants. | The discretization of segments is limited and cannot be handled irregularly. It also requires a lot of experimental field data. | [27] |
| MIKE 11 | Applied in the evaluation of water quality in water bodies affected by dumping and eutrophication. | Several variables associated with heavy metals are not available. Detailed and not available data about the water body are required. | [28] |
| EFDC | 3-D simulation of hydraulic data from different water bodies. EFDC can be coupled with other model systems. | Simulation only admits hydrodynamic transport for pollutants. It lacks graph interphase as output. | [29] |
| WASP | 3-D simulation for prediction of affectation in the water quality for any hydric system. Previous studies show its application in heavy metals. It offers 2-D graph outputs. | It requires external sources or previous results for hydrodynamic calculus. | [30] |
[i]Input data from model [31]