A coupled model of hydrodynamics circulation and water quality applied to the Rio Verde reservoir , Brazil

This work applied the coupled horizontal two-dimensional hydrodynamic circulation model (2DH) and the vertically integrated water quality model for non-conservative and passive scalars to the Rio Verde reservoir in the state of Paraná, Brazil, to simulate flow, temperature and water quality parameters. The water quality model uses the same spatial grid applied for the hydrodynamics model. Flow velocities and turbulence coefficients previously defined in the hydrodynamics model can be used directly in the water quality model. Modeling results were compared to field data for a period of 308 days, from February 27 2010 and December 31 2010. Field data included water and air temperature, wind, relative humidity, radiation, discharges and concentration values of some substances in the tributaries of the reservoir. The results indicated that hydrodynamic circulation and, consequently, horizontal transport, are strongly dependent on the wind. Inflows/outflows generated a localized circulation. The results for the temperature and dissolved oxygen concentration were compared with field measures and a satisfactory consistency was achieved. Despite the errors associated with boundary conditions, the models demonstrated their potential to adequately simulate the data set collected from the reservoir.


INTRODUCTION
Reservoirs are complex environments which create interdependent relationships between their tributaries (rivers and streams that contribute for reservoir discharge) and the watershed environmental conditions, which are strongly influenced by water resources uses.By reducing rivers' velocities with dams, lentic environments can be created which can in turn drastically alter the original conditions of the rivers with greater flow (lotic environments).In this context, the importance of limnological studies aiming at understanding the transport and mixing mechanisms of substances in lentic environments is emphasized.
Quantitative water quality modeling is a major challenge of these studies, due to the need for interdisciplinary approaches that involve advective-diffusive mechanisms and physical, chemical and biological processes in substance transportation.As a consequence, the development of mathematical models to quantify these processes has been the subject of intense research in the last decades.System complexity should be taken into account during the modeling process, ensuring that it has the capacity to model the environmental behavior and evaluate its several effects.Numerical models have been used as a scientific and managerial tool for analyzing temporal-spatial distributions of non-conservative water quality parameters.Among many works, the following deserve to be mentioned: Deus et al. (2013) show the application of the CE-QUAL-W2 model in the Tucurui reservoir in Brazil for different pisciculture scenarios.Using the same model, Kuo et al. (2006) simulated water quality changes at temperate and sub-tropical reservoirs in Taiwan.Liu et al. (2008) published a study describing a 2-D coupled model of hydrodynamics and water quality applied to Yuqiao reservoir in China, developed to compute the hydrodynamic field and the variations of total phosphorus and total nitrogen.
The Rio Verde reservoir is located in the Metropolitan Region of Curitiba (RMC), in the state of Paraná, Brazil.Based on collected data included in the work by Carneiro et al. (2014), one can notice that the Rio Verde reservoir presented a moderately degraded environment.In 2005, an algae bloom of Cylindrospermopsis reciborskii was recorded, reaching 96,489 cells.mL - .The Rio Verde reservoir has a dynamic system of stratification characterized by a reasonable stratification in the summer and a mixed water column in the winter.Due to the shallow depths, it is expected that the currents flowing can be well represented by depthaveraged variables.
Therefore, the purpose of this study was to use the hydrodynamic and water quality models of the SisBAHIA ® , in coupled form, to simulate flow, temperature and water quality parameters in the Rio Verde reservoir.The advantage of coupling the two models, as developed in the current work, appears in the determination of the flow velocities and turbulence coefficients, Rev. Ambient.Água vol.13 n.6, e2244 -Taubaté 2018 which is done previously in the hydrodynamic model, and can be used directly in the water quality model.For this reason, the size of the system to be solved and, consequently, the computational effort, are reduced when compared to other schemes.
The SisBAHIA ® (in Portuguese, Sistema Base de Hidrodinâmica Ambiental) is the Hydrodynamic Environmental System developed by the Coastal and Oceanographic Engineering Department, Oceanic Engineering Program, Federal University of Rio de Janeiro (COPPE/UFRJ).In the present work, only the water quality model of SisBAHIA ® is shown.Further information about the hydrodynamic model of SisBaHiA ® can be obtained through: www.sisbahia.coppe.ufrj.br.

MATERIALS AND METHODS
The water quality model of SisBaHiA ® (MQA) is an Eulerian advective-diffusive transport model vertically integrated for non-conservative scales; in other words, the variability of substances concentrations is a result of physical, chemical and biological processes.
In the discretization of time and space, the SisBaHiA ® adopts the finite difference method and the finite element method, respectively.The MQA uses the same spatial grid applied for the hydrodynamics model and different time-step lengths can be employed in analyses.Flow velocities and turbulence coefficients, previously defined in the hydrodynamics model, can be used directly in the water quality model.The advective mechanism is related to flow velocities and can be estimated through hydrodynamic models, which represent an important previous procedure for water quality models' input data.The time step of the two-dimensional hydrodynamic model was 30s with an average Courant number of 4.0.The time step of the MQA was 60s.
The MQA is able to simulate the oxygen, nitrogen and phosphorus cycles, as well as temperature and salinity parameters due to their strong influence in kinetic processes (Sellers, 1965).Overall, 11 water quality parameters can be evaluated, namely: Salinity, Temperature, Dissolved Oxygen (DO), Biochemical Oxygen Demand (BOD), Organic Nitrogen, Ammonia Nitrogen, Nitrate Nitrogen, Chlorophyll_a, Herbivore Zooplankton, Organic Phosphorus and Inorganic Phosphorus.The MQA uses an Eulerian approach by the finite element mesh, previously defined in the hydrodynamic model.With this condition it is possible to share the same flow velocity components, wind velocity and geometric characteristics.The spatial discretization employs, mainly, squared biquadratic finite elements (Rosman, 2016).

Mathematical model
In the MQA, the mass-balance equation for a non-conservative substance is applied considering depth-integrated Cartesian System (x, y) aligning in the east, north and vertical directions for three transport mechanisms: the advective term, the diffusive term, and the kinetic processes.The mass-balance equation is expressed by Equation 1, as defined by Sheng and Villaret (1989): In Equation 1, Cm is the concentration of m substance (mg L -1 ); t is the time (s); Ui represents the depth-averaged components of the horizontal velocity (U and V) (m s -1 ); H is the water depth (m); Dij is the turbulent viscosity coefficient of mass (m 2 s -1 ); δjk is the Kronecker delta; Λk represents the widths of the spatial and temporal Gaussian filters; and Rm represents the kinetic processes of substance production and consumption.In Equation 1, i, j = 1,2 and k = 1, 2, 3, with k = 3 corresponding to time t.The following interpretation is valid for the index coefficient Cm: C1 = ammonia nitrogen (mg N L -1 ), C2 = nitrate nitrogen (mg N L -1 ), C3 = inorganic phosphorus (mg P L -1 ), C4 = Herbivore Zooplankton (mg C L -1 ), C5 = biochemical oxygen demand (mg O2 L -1 ), C6 = dissolved oxygen (mg O2 L -1 ), C7 = organic nitrogen (mg N L -1 ), C8 = organic phosphorus (mg P L -1 ), C9 = Chlorophyll_a (μg L -1 ), CT = temperature ( 0 C) and CS = salinity (psu).
The MQA is coupled with the SisBaHiA ® hydrodynamic model to provide the advective and diffusive terms in the mass-balance equation.The kinetic processes in MQA were obtained according to the Equations 2, 3, 4, 5, 6, 7, 8, 9 and 10 below (Sheng et al., 1996).C1: Ammonia nitrogen C5: BOD -Biochemical oxygen demand The numerical implementation of the two models, hydrodynamic and water quality, are not discussed here; the interested reader is referred to Rosman (2016) for additional information concerning the hydrodynamic model.The numerical model developed for the advective and diffusive transport is described in detail by Cunha et al. (2002).

Description of the Rio Verde Reservoir
The Rio Verde reservoir is located in the Metropolitan Region of Curitiba (RMC) and inserted in the Rio Verde watershed, which has a total area of 242 km 2 .The reservoir has a plan area of 5.971 Km 2 , an average volume of 25.6 hm³, an average water depth about 5.6 m and receives the water contribution of 14 sub-watersheds.The Rio Verde reservoir is oriented in a Northeast (NE)-Southwest (SW) direction, with a length of approximately 7,500 m and a maximum width of 1,300 m.The ratio between the depth and width of the reservoir is about 1:230 and between the depth and length is about 1:1330; these features suggest a predominance of horizontal velocities.
Figure 1 shows the bathymetry data and also the water quality and discharge points used for measurement.The main tributary of the reservoir is the Verde River (F4), accounting for about 70% of the total discharge flow; the other 13 smaller tributaries contribute with low flow levels.Figure 2 shows the daily discharge of each tributary from March 2010 to December 2010, estimated by the SWAT model (Soil and Water Assessment Tool).In relation to outflows, two points can be set: one close to the dam, which represents a REPAR (President Getúlio Vargas Refinery) license, with value of 0.83 m 3 /s; and another at the spillway, which represents difference between the sum of inflows and outflows (Carneiro et al., 2014).Meteorological data were obtained from a station installed near the dam (Latitude 25º31´36,83" S and Longitude 49º31´39,07" W).This station provided wind direction, wind velocity, air temperature, solar radiation, relative humidity and precipitation from March 2010 to December 2010.The daily maxima values of solar radiation showed a gradual decrease from March to July, reaching 457 W m -2 in May (Figure 2); from August to December the values Rev. Ambient.Água vol.13 n.6, e2244 -Taubaté 2018 increased, reaching 851 W m -2 in November.In Figure 2, the air temperatures data show a variability related to seasons, with minimum moving average of 12.6°C in August and maximum moving average of 20.7°C in April.With respect to the relative humidity, Figure 2 shows a small moving average variation over 2010, with values between 74.1% and 88.2%.Carneiro et al. (2014) show that only the Verde River (F4) and two tributaries (TE10 and TD4) account for 90% of the pollutant load that arrives at the Rio Verde reservoir.Figure 3 shows the concentration values of eight parameters in three tributaries (F4, TE10 and TD4) from March 2010 to December 2010: water temperature, ammonia, nitrate, organic nitrogen, organic phosphorus, inorganic phosphorus, DO and BOD.
The Rio Verde reservoir presents a dynamic system of stratification characterized by a reasonable stratification in the summer and a mixed water column in the winter.Figure 4 shows the temperature profiles measured at R1 and R4 stations (see Figure 1) from March 2010 to December 2010.At the R4 station between April and September, the water column shows no stratification and can be considered well mixed; from October the water column begins to stratify with the increasing air temperature, reaching largest gradients in December.The station R1 is located in the shallow portion of the reservoir, and it is not possible to observe any stratification phenomena there.
With respect to the wind characteristics, the measured data between March and December 2010 demonstrate a predominance of NE direction, which coincides with the alignment of the reservoir axis.The wind flow represents the major influential variable for hydrodynamic circulation and mass transport in this reservoir.In two areas the hydrodynamic flow and mass transport are dominated by inflows and outflows: one near the dam, where the main water intake and the spillway are located; and the other near the Verde River, the main tributary.Outside these areas, circulation and transport are determined by the wind.

RESULTS AND DISCUSSION
Two-dimensional hydrodynamic models of Rio Verde reservoir were developed using SisBaHiA ® .The spatial discretization in the horizontal x-y plane was carried out through subparametric Lagrangian quadrilateral finite elements with 9 nodes.In a sub-parametric element, the geometry of the element is linear and is defined only by the vertices.However, the variables are quadratic, and in addition to the values at the vertices, the middle values at each side are also needed and, in the case of quadrilaterals, a node at the center of the element is also necessary.
The mesh consists of 507 elements and 2402 nodes.The simulation was performed for a period of 308 days, from February 27 th 2010 to December 31 th , 2010, with a 60s time step.These dates were chosen due to available measurements of water quality parameters, made at four locations in the reservoir, which enabled the calibration process in MQA evaluation.
In the numerical model, wind conditions are considered unsteady and spatially homogeneous; the input data used in the model were the time series of wind speeds and directions measured at the meteorological station in the reservoir.At the open boundaries domain, the elevations were imposed; at the land boundary, except for the nodes corresponding to the rivers, all nodes were considered impermeable with null tangential and normal velocity components imposed.
The bottom friction coefficient can be written in terms of the Chezy coefficient, which depends on the amplitude of the equivalent bottom roughness.The amplitude of the equivalent bottom roughness, ε, was defined according to the bottom sediment characterization and Rev. Ambient.Água vol.13 n.6, e2244 -Taubaté 2018 distribution (Rosman, 2016).Silty-clay sediments are predominant in the reservoir (ε ≈ 0.015m) and in regions with large amount of submerged trees, bottom roughness amplitude values were increased (with values around 0.130 m).
Figure 5 shows the hydrodynamic model results for two situations: in June 26 th , 2010, with an average wind speed of 1,81 m s -1 , an average wind direction of 73°, and a main river discharge of 2,24 m 3 /s; and in October 18 th , 2010, with an average wind speed of 4,84 m s -1 , an average wind direction of 22°, and a main river discharge of 1,63 m 3 s -1 .The velocity fields show a strong relationship with the local wind pattern through the presence of a clockwise vortex formed in the region near the dam, which intensifies depending on the wind pattern.On June 26 th , 2010, there was a northeasterly wind at an average of 1.81 m s -1 .The currents follow the wind direction, forming a large clockwise vortex in the region near the dam.In the regions near the spillway and the Verde river, there is an intensification of the currents due to local effects of inflows/outflows.On October 18 th , 2010, the wind speed is stronger, at an average of 4.84 m s -1 .Under these conditions, two vortices were formed: a large vortex, clockwise in the central region of the reservoir, and a smaller one, also clockwise.There was an intensification of the currents in the region near the Verde river, caused by the flooding observed during this period.Taking into account these patterns, it can be stated that the hydrodynamic circulation, and consequently the horizontal transport, is strongly dependent on the wind and that inflows/outflows generate a localized circulation.From March 2010 to December 2010, data related to Temperature and DO were collected.The measurements were made at R4 station in reservoir, see Figure 1; the measurements were obtained at a several positions along the water depth.The analyses were performed for the same dates as that of the hydrodynamic simulation, between February 27 th 2010 and December 31 th , 2010.This article shows just the results for the Temperature and DO, which have measured data for comparison.
The turbulent diffusion coefficient used from the water quality model are: Dxx (m 2 s -1 ) = 2.0, Dyy (m 2 s -1 ) = 0.5, with a 60s time step.Table 1 lists the model coefficients and constants employed in the numerical simulation.The values of the parameters and constants are welldefined in the literature and were used in the application considering that there are no available studies for water quality modeling in Rio Verde reservoir.This is the main difficulty in performing the analysis proposed in this article.The field data are not sufficient for a complete calibration, and some adjustments were carried out manually to obtain numerical results close to the measured data, within a specified variation limit; some adjustments proved to be necessary, mainly in the reaeration and deoxygenation coefficients.In order to quantify the accuracy of fit, correlation coefficient R 2 of the regression line between the model results and the field data at R4 station was computed.The best calibration required a R 2 value as close to 1.0 as possible.
The temperature is used as the calibration parameter for variables related to advective and diffusive transport.The long-term temperature data simulated with SisBaHiA ® were compared with average values of vertical profiles, measured at R4 station from March 2010 to December 2010.The model reproduces the average temperature variation and, according to Figure 6, it is possible to state that variable boundary conditions are adequately reproduced.The R4 is not under the direct influence of the rivers; consequently, there is a better agreement between the results.Figure 6 shows the DO concentrations measured at the R4 station and computed numerically by SisBAHIA ® .The differences between these results showed that SisBAHIA ® can predict adequately the variations of DO in the process of oxidation of organic matter present in the reservoir.The average relative errors between the model results and the field data at R4 station are: 0.095 for temperature and 0.406 for DO. Figure 6 also depicts the regression line; the correlation coefficient values (R 2 ) were 0.8973 for the temperature, and 0.4575 for the DO.This indicated good agreement between simulated and field data for the temperature.However, DO concentrations do not indicate a good fit.For a better calibration of the model, a larger temporal series would be necessary.

CONCLUSIONS
The results presented in this paper show that the environmental modeling of the Rio Verde reservoir is appropriate and can be used quantitatively to study the phenomenon of eutrophication.Whereas the advective and diffusive transport is independent from the parameter modeled, the correct definition of the chemical, physical and biological processes involved in each water quality parameter, made it possible to characterize the distribution of these parameters in the Rio Verde reservoir.Based on the results of the hydrodynamic model, it was possible to note the significant influence of the wind in the hydrodynamic circulation.Thus, the importance of the proper use of wind data when modeling reservoirs becomes clear.Inflows/outflows generate impacts in circulation only in their immediate surroundings and do not contribute to the circulation of the reservoir as a whole.The results for the temperature and kg -Phytoplankton growth rate (d -1 ) 2.0 ksN -Half-Saturation Constant for Nitrogen (µgN L -1 ) 25.0 ksP -Half-Saturation Constant for Phosphorus (µgP L -1 ) 20.0 Is -Optimal light level (ly d -1 ) 250.0 ksa -Half saturation constant for predation of zooplankton on algae (µgChla L -1 ) 15.0 kgz -Grazing rate (m 3 gC -1 d -1 ) 1.0 kra -Phytoplankton losses due to respiration and excretion (d -1 ) 0.05 krz -Zooplankton losses due to respiration and excretion (d -1 ) 0.01 Er -Efficiency predation of zooplankton on algae 0.6 kgzc -Predation losses rate (d -1 ) 0.01 f0n -Fraction of death and respiration of phytoplankton recycled to organic nitrogen 0.5 Kam -Half-saturation constant for preference of Ammonia Nitrogen (μgN L -1 ) 50.0 roc -Oxygen amount consumed in the decomposition of one gram of organic carbon (gO gC -1 ) 2.67 rpa -Ratio of Phosphorus to chlorophyll in the phytoplankton (gP gChla -1 ) 1.0 rna -Ratio of Nitrogen to chlorophyll in the phytoplankton (gN gChla -1 ) 7.2 rca -Ratio carbon/chlorophyll in phytoplankton cells (gC gChla -1 ) 50.0 fD5 -BOD fraction dissolved in the water column 0.5 fD7 -Fraction of dissolved organic nitrogen in the water column 1.0 fD8 -Fraction of dissolved organic phosphorus in the water column 0.85 fop -Fraction of dead and respired of phytoplankton in phosphorus cycle 0.5 k12 -Nitrification coefficient (d -1 ) 0.1 k2D -Denitrification coefficient (d -1 ) 0.1 k71 -Organic nitrogen mineralization coefficient (d -1 ) 0.03 k83 -Organic phosphorus mineralization rate (d -1 ) 0.03 ka -Reaeration coefficient (d -1 ) 1.38 kD -Deoxygenation coefficient (d -1 ) 0.2 kDBO -Half saturation constant for oxidation of BOD (mgO2 L -1 ) 0.5 knit -Half-saturation constant for DO limitation in the nitrification process (mgO2 L -1 ) 0.5 kno3 -Half-saturation constant for DO limitation in the denitrification process (mgO2 L -1 ) 0.1 kea -Phytoplankton mortality rate (d -1 ) 0.1 kez -Zooplankton mortality rate (d -1 ) 0.005 SOD -Sediment oxygen demand (gO2 m -2 d -1 ) 1.0 vs3 -Organic matter settling velocity (m d -1 ) 1.0 vs4 -Phytoplankton settling velocity (m d -1 ) 0.05 vfr -Inorganic sediment settling velocity (m d -1 ) 0.02

Figure 1 .
Figure 1.Map of Rio Verde reservoir bathymetry, water quality measure points and discharge measure points during 2010.

Figure 2 .
Figure 2. Tributaries daily discharge (m 3 s -1 ), daily average of solar radiation (W m -2 ), air temperature and relative humidity measured at the meteorological station of the Rio Verde reservoir from March 2010 to December 2010.The continuous line represents the moving average for 30 days.

Figure 3 .
Figure 3. Concentrations values of some substances in the tributaries of the Rio Verde reservoir.

Figure 6 .
Figure 6.Temperature and DO concentration measured at the R4 station and computed numerically by SisBAHIA ® .
(1 −   )   9 +  oc  ca    9 +      4  6 =   (  −  6 ) −    5 −    12  1 +  oc  ca    9 −        9 −      4 −   (1 −   )   9  4 +        9 + ( The coefficients used in the model are given in Table 1.In the MQA two different types of horizontal boundaries are considered: land boundaries and open boundaries.The land boundaries, in general, represent the water body margins and points with inflows or outflows, such as rivers, streams, surface water abstraction points, among others.The open boundaries represent the water domain limits of the study area.The prescription of normal fluxes in this model is associated with land boundaries.Along open boundaries, it is usual to neglect the diffusive fluxes; the model computes the mass balance equation with no diffusive terms. 9 =    9 −    9 −    9 −    9 −     9 (10) Rev. Ambient.Água vol.13 n.6, e2244 -Taubaté 2018

Table 1 .
Parameters employed in the water quality model and values adopted in the Rio Verde reservoir analysis.