Comparison of 1D and 3D reservoir heat transport models and temperature effects on mass transport

One and three-dimensional heat transport models are compared in a dendritic reservoir in Brazil. We estimate the periods of temperature stratification for both models using physical indices and temperature gradients. The three-dimensional model reproduces more accurately the water column temperature profiles, however with focus on the physical indices (Wedderburn Number and Lake Number) similar results were obtained with both models. Secondly, we investigated the effects of temperature stratification on substance mass transport using the three-dimensional model. The advective and dispersive transport for a tracer released in a river and in a side arm of the reservoir were quantified. We identified that considering the effects of temperature stratification increased the horizontal advective transport up to a maximum of 40% increase for the tracer released in the river, and 9% for the side arm. In relation to dispersive transport there was a decrease in transport due to temperature stratification, and no pattern was discernible for the side arm tracer modeling.


INTRODUCTION
Lakes and reservoirs vary widely in their size and number. They comprise approximately 0.02% of the water in the hydrosphere (MARTIN; MCCUTCHEON, 1999), in contrast of 0.00009% in rivers (WETZEL, 1983). Usually, it is distinguished between reservoirs and lakes. Latter is considered natural and the former is typically a result of engineering projects (MARTIN; MCCUTCHEON, 1999) with many purposes, such as water supply, irrigation, flood control, navigation and power generation (UNESCO, 2009) and often used for multiple purposes (MARTIN; MCCUTCHEON, 1999).
Construction of large reservoirs in Brazil began in the early 1900s (SOARES et al., 2008). Brazil is a reservoir-oriented country (TUNDISI et al., 1993). According to the National Energy Balance of 2013, hydropower energy corresponds to 70.1% of the Brazilian demand. Between 2011 and 2012 there was an increase of 2.2% in hydroelectric capacity (BRASIL, 2013). Reservoirs used for hydropower generation may have multiple purposes, as water supply, irrigation, and recreation (SOARES et al., 2008), therefore both, water quantity and quality are important for their operation.
Besides hydropower reservoirs, the number of drinking water reservoirs is also increasing. According to ANA (ANA, 2013) -the Brazilian Water Agency, flow consumed by multiple purposes in 2010 had the proportion of: 72% to irrigation, 11% to animal consumption, 9% to urban consumption, 7% to industry uses and 1% human rural supply. Surface waters, which include rivers, lakes, and reservoirs, are still the most common used sources for drinking water supply -56% of the municipalities use surface waters as their main water source (GALLI; ABE, 2010).
Construction of reservoirs may impact physical, chemical and biological characteristics of the water body by affecting residence time, water temperature, stratification, reduction in mixing characteristics (less turbulent waters) and often a decrease in particles and turbidity (FRIEDL; WUEST, 2002). Some water quality issues due to river damming are eutrophication (DOAN et al., 2015), dissolved oxygen depletion (SCHLADOW; HAMILTON, 1997;GIANNIOU, 2003;BELL et al., 2006), manganese release (BERTONE et al., 2015) and greenhouse gas emission to the atmosphere (DEEMER et al., 2016). Knowledge of the physical processes that occur inside the reservoir and its driving forces are important for the management of these systems (CURTARELLI et al., 2014). From the environmental perspective, studies about the impacts of dams are required for the planning and licensing of these engineering structures. Already built reservoirs also require knowledge of the physical-chemical processes in relation to operation and maintenance.
Stratification controls the distribution of heat and substances in lakes and reservoirs (ESTEVES, 1989). Therefore, water column temperature is an important parameter in studies related to circulation and water quality in reservoirs. In the last years, attention has been given also to the horizontal distribution of temperature, that may also affect transport, and act as a source of substances to deep waters. This is due to horizontal processes that may induce circulation between shallow and deep water. This circulation could be a source of dissolved oxygen (HOHMANN et al., 1997), phosphorus (JAMES; BARKO, 1991), for example, which could interfere in the structure of the water column and water quality parameters.
Heat transport in lakes may be modeled by one, two or three-dimensional models. One-dimensional models have been widely used to simulate seasonal variations in reservoirs, for the identification of mixing periods and sometimes coupled to water quality modules (MARTIN; MCCUTCHEON, 1999;WEINBERGER;VETTER, 2012). These models assume horizontal homogeneity, and vertical turbulent diffusivity being modeled as a relation of wind speed at the surface and usually not being uniform (JASSBY;POWELL, 1975;SUNDARAM;REHM, 1971;HENDERSON-SELLERS, 1986). Most models consider a no flux condition for heat energy at the bottom (GIANNIOU; ANTONOPOULOS, 2007;BABAJIMOPOULOS;PAPADOPOULOS, 1986), although that could not be true for shallow reservoirs (WETZEL, 1983). At the surface the energy flux is parametrized (EDINGER et al., 1968;SUNDARAM;REHM, 1971) or calculated term-by-term (HENDERSON-SELLERS, 1986). Some of the most known models are: DYRESM (WEINBERGER; VETTER, 2012; HOCKING; PATTERSON, 1991;IMERITO, 2007), SIMSTRAT GOUDSMIT et al., 2002), FLake (MIRINOV et al., 2003) and GLM (HIPSEY et al., 2014). Polli and Bleninger (2018) recently developed a one-dimensional model, MTCR-1, showing that errors compared to field measurements were comparable to some of these well-known models.
Laterally averaged two-dimensional models simulate hydrodynamics, mixing, and water quality changes in the longitudinal and vertical direction (MARTIN; MCCUTCHEON, 1999). They are applicable to stratified reservoirs where lateral variations are small (MARTIN; MCCUTCHEON, 1999). Two-dimensional models solve the two-dimensional laterally averaged equations of continuity, momentum and transport. CE-QUAL-W2 is one well known laterally averaged model and has been widely used in problems related to thermal stratification (MODIRI-GHAREHVERAN et al., 2014;ZOUABI-ALOUI;GUEDDARI, 2014;BONALUMI et al., 2012).
Three-dimensional models solve the full continuity, momentum and transport equations. Some well-known three-dimensional models are ELCOM and Delft3D. Both models solve the Reynolds Averaged Navier Stokes equation, assuming hydrostatic conditions in the vertical and applying the Boussinesq approximation (DELTARES, 2014;LEE et al., 2013). The numerical approach used in ELCOM is the finite volume method (CHUNG et al., 2009) while in Delft3D a finite difference method is used (DELTARES, 2014). Both models use a staggered grid (DELTARES, 2014;LEE et al., 2013) which defines velocities on cell faces and scalars on cell centers. Both models follow the stability criterion of Courant-Friedrichs-Lewy (DELTARES, 2014;LEE et al., 2013). Delft3D has five different heat flux models (DELTARES, 2014), while ELCOM has one heat flux model that calculates all heat fluxes at the air-water interface separately (LEE et al., 2013). Delft3D has four models for the turbulence closure for vertical viscosity and diffusivity (from simple uniform value type models up to a k-ε model that solves two transport equations for production and dissipation Polli; Bleninger 3/19 of turbulent kinetic energy (DELTARES, 2014). ELCOM uses an approach derived from the mixing energy budget used in 1D lake modeling (LEE et al., 2013). ELCOM estimates an energy budget between the energies available for mixing, required for mixing, and dissipated energies, based on the 3D turbulent kinetic energy transport equation (HODGES et al., 2000).
Considering the huge number of reservoirs, and the large variety of available modeling tools, questions arise on which model to use for which purpose, and optimal model setups, considering the available heat-transfer formulations, as well as relevant hydrodynamic features, originating from stratified flows. Some relevant questions for the management of reservoirs are: Does the reservoir stratify? If yes, when and for how long? What are the consequences of stratification? Does it affect water quality? To answer these questions, usually numerical models are applied, since they complement the understanding of physical processes (MARTIN; MCCUTCHEON, 1999) and allow simulations of real and hypothetical scenarios where various processes and forcings can be applied simultaneously or separately (CURTARELLI et al., 2014). In general, there is a lack of measurements for Brazilian reservoirs, therefore, the implementation of numerical models is complex. On the other hand, model applications showed to be very useful to understand the dynamics of the system. One-dimensional temperature profile models allow simulations of long period, with high temporal and spatial resolution, and low processing time -but they do not include horizontal advection, nor allow for analysis of different locations in the reservoir. Three-dimensional models are usually adopted for short time periods and with a coarse grid, due to their elevated processing time and post-processing needs. Seasonal variations are thus not always completely analyzed and, from the operating or management point of view they are not practical and efficient, due to complexity and long simulation times. But 3D models can simulate the entire system and provide local analysis.
In this paper, it is analyzed heat transport due to its role on the dynamics of reservoirs and implications for hydrodynamics and water quality. We compare one and three-dimensional heat transport models for Vossoroca reservoir. Then, it is investigated and quantified (advective and dispersive transport) if temperature may enhance or be responsible for mass transport considering a three-dimensional simulation of a tracer continuously released in a river and in a side arm of that dentritic reservoir. This paper describes mainly for two contributions related to the study of heat and mass transport in reservoirs. The first one, to answers the often arising question "when do 1D models are sufficient to describe reservoir dynamics", even though knowing that the processes in reservoirs are mainly three-dimensional. We showed that for heat transport modeling a in each case a one-dimensional model describes the thermal structures (i.e. stratification) well enough. The second contribution answers the question "what is the role of reservoir stratification on mass transport", and "how the system responds when temperature is considered and not considered in the modeling". We showed that indeed temperature stratification plays a significant role for the dynamics of the horizontal mass transport, and its vertical position (layering) in the reservoir.

Site description
Vossoroca reservoir ( Figure 1) is a monomictic reservoir, located in Tijucas do Sul, near to the city Curitiba, in Brazil (49ºW 25º52'S) at 800 m above sea level. It was formed in 1949 by the impoundment of São João river, and its function is to regulate the flow to Chaminé hydropower plant (with capacity of 18 MW). The surface area is 3.3 km 2 and the volume 35.7 10 6 m 3 . It's maximum and mean depths are 17 and 8 m, respectively and the residence time, 117 days. The rivers flowing to Vossoroca reservoir (see Figure 1) are: São João river (56% of the total flow), São Joãozinho river (28% of the total flow) and Vossoroca river (16% of the total flow) (BERNARDO, 2013). The climate of this region is humid subtropical and its annual average temperatures varies from 18 °C and 22 °C. Vossoroca reservoir is monitored since 2012. Its monitoring includes water temperature, meteorological measurements and inflows/outflow. The next section will present measured data in the reservoir and required data to run the models. Figure 1 shows the location of the monitoring stations. For this study, the period of measurements used is from July to October (07/01/2012 to 10/02/2012), that represents an initially well-mixed period, followed by a period of strong stratification and continuous water level variation.

Field measurements
Temperature measurements were taken at the deepest point of the reservoir using 7 temperature sensors. The depths measured (from surface) are: 1, 3, 5, 7, 9, 11 m and 1 m measured from the bottom, with 15 minutes time step (MÄNNICH, 2013). Recently, this configuration has been changed and 30 temperature sensors were deployed at this point. Figure 2 shows temperature profile measurements at the floating platform location. During Winter in the Southern Hemisphere the reservoir is well mixed, with water temperatures around 14ºC. Subsequently, the reservoir warms at the surface, developing a thermal stratification after August.
Meteorological data is measured with 2 min time step at an onshore weather station. The meteorological weather station measures precipitation, air temperature, relative humidity, solar radiation, pressure, wind speed and wind direction (MÄNNICH, 2013). Figure 3 shows the meteorological information required to run both one and three-dimensional models.
Flow and water level in the reservoir are monitored by COPEL (electrical company of Paraná State). For the period of interest of this work, the mean inflow was 2.5 m 3 s -1 , and outflow 3.7 m 3 s -1 , consequently resulting in a decrease of the water level. Figure 4 shows the discharge and inflow temperature for the period of interest. Usually, inflow temperature is lower than the water column temperature in the reservoir -this way, most of the time, the inflows enter the reservoir as an underflow, which was also verified by CTD measurements in the reservoir.

Model description
One-dimensional model: MTCR-1 MTCR-1 model solves the 1D heat transport equation (POLLI; BLENINGER, 2018; POLLI; BLENINGER, 2015). One-dimensional heat transport models assume the vertical direction being representative of the heat transport, therefore, the main assumption is the homogeneity in the horizontal direction. Heat transfer across the air-water interface is calculated according to equilibrium temperature (EDINGER et al., 1968) and a no heat flux condition at the bottom layer. The eddy diffusivity coefficient is calculated according to Sundaram and Rehm (1971), and it is a function of the Richardson number. The model simulates inflows and outflows as a source of heat and turbulence to the reservoir, and there is a verification of instabilities of the temperature profile at each time step. Also, it is possible to simulate high temporal and spatial resolution in MTCR-1 model, with low processing   time. The heat transport equation is solved with an implicit finite volume method and programmed in Fortran 95. MTCR-1 was already calibrated for Vossoroca reservoir (POLLI; BLENINGER, 2018; POLLI; BLENINGER, 2015) and for Harp Lake, in Canada (GUSEVA et al., 2017). More information about the model may be found in those publications.
The grid in the MTCR-1 model was defined with a maximum of 85 fixed layers, which correspond to 0.2 m thickness of each layer for Vossoroca reservoir. The model considers the variation of the reservoir area with depth. The UNESCO equation is used to calculate water density for fresh water. Solar radiation is absorbed as a function of depth -Secchi depth is a calibration parameter for the model. Secchi depth was set to 1.8 m, which is close to the measurements at the deepest point of the reservoir. The mean absolute error (MAE) and standard deviation (SD) were calculated for estimation of errors in relation to the measurements. The time step for the simulation was set to 2 min. In July, Vossoroca reservoir was mixed, therefore the initial condition of the temperature set to the measured temperature of 14.5 °C. Three-dimensional model: Delft3D Delft3D (DELTARES, 2014) is an open source model and has been widely used to solve problems related to reservoirs and thermal stratification (SOULIGNAC et al., 2017;WAHL;PEETERS, 2014, SMITS et al., 2009. The hydrodynamic module Delft3D-FLOW simulates the two-dimensional (depth averaged) and three-dimensional unsteady flow and substance transport resulting from tidal and/or meteorological forcings, and includes the effect of density differences (due to non-uniform temperature and salinity). The hydrodynamic modeling system solves the unsteady shallow water equations for an incompressible fluid. The system of equations is composed of the equations of motion, the continuity equation and the transport equation for conservative constituents (DELTARES, 2014).
The horizontal grid was defined by 342 points in M-direction (x direction in the Cartesian plan) and 347 in N-direction (y direction in the Cartesian plan) -the size of the Cartesian grid cells are 15 m x 15 m. In the vertical direction, a Z-model was used with 20 layers, since this configuration is recommended for a better representation of thermal stratification in lakes and reservoirs. A σ layer model was tested, but indeed resulted in increased mixing in the water column. The UNESCO equation was used to calculate water density and salinity was set to zero in the model.
For the turbulence closure scheme, a k-ε model was used to calculate the vertical eddy viscosity, and vertical eddy diffusivity. For the background vertical viscosity and diffusivity, the coefficients were considered as calibration coefficients, and fixed during the simulation. The background horizontal viscosity and diffusivity were also used as calibration coefficients and fixed during the simulation. Table 1 shows the parameters used for turbulence parameter calibration.
Delft3D considers that there is no heat flux between the bottom of the reservoir and the water. At the surface, five heat flux models are available to estimate air-water fluxes. We compared Murakami and Ocean models to define the most adequate model to be used for Vossoroca reservoir. Both models need Secchi depth for solar radiation absorption, used as a calibration parameter. Dalton number is used for calibration of the evaporative heat flux in both models. In the Ocean model, Stanton Number for heat convection is a calibration parameter. Table 1 shows the parameters calibrated for surface heat flux for Vossoroca reservoir modeling. The calibration parameters were defined as the ones that best fitted measured temperature profiles at the floating platform location. The criterion used was the smallest mean absolute error (MAE) together with standard deviation (SD).
The time step of the simulation was set to 3 s, respecting Courant-Friedrichs-Lewy number (CFL). Initial conditions were defined according to the measurements and considered as uniformly distributed values over the whole domain. In July, Vossoroca reservoir was well mixed and the initial temperature was 14.5 °C. During the period of interest, the water level varied in Vossoroca reservoir around 2 m -the initial condition was set to -0.42 m.

Physical Indices
The modeling results (i.e. temperature profiles) of both models were analyzed and compared quantitatively, using physical indices to classify the mixing characteristics. Lakes and reservoirs may be classified according to its periods of mixing and stratification. Read et al. (2011) developed the software Lake Analyzer, that quantifies important parameters related to these periods in lakes and reservoirs, based on measurements. The characterization of the system according to mixing and stratification periods can be performed using Lake Number and Wedderburn Number. The Wedderburn Number is a relationship between buoyancy and mixing induced by wind shear, calculated by (IMBERGER; PATTERSON, 1989): in which g' is the reduced gravity (g'=gΔρ/ρ 0 ), g is the acceleration of gravity, Δρ is the density difference between the hypolimnion (ρ h ) and epilimnion (ρ e ), z e is the epilimnion thickness (m), L is the reservoir fetch length (m). If W<<1, the increase of the mixed layer depth is dominated by the internal production of turbulence. In this case, upwelling can occur upwind and strong gradients at downwind (IMBERGER; PATTERSON, 1989). On the other  , 1989). Stratification is strong and the mixed layer will deepen slowly (READ et al., 2011). The Lake Number describes relevant processes for mixing in lakes induced by wind action, according to (READ et al., 2011) ( ) where S t is the schmidt stability that relates the resistance to mechanical mixing due to the potential energy on the stratification of the water column, z h is the Hypolimnion thickness (m), A s is the surface area of the reservoir (m 2 ) and (Read et al., 2011). If L N >1, stratification is strong and any disturbance produced by the wind at the surface is minimized. Otherwise, if L N <1 indicates weak stratification and potential to mixing (IMBERGER; PATTERSON, 1989).

RESULTS ANS DISCUSSION
The results are presented in the following three sections. In the first section, two heat flux models available in Delt3D are compared for the identification of the most adequate model to use for Vossoroca modeling. The second section presents the one and three-dimensional model results. The third section shows the effect of temperature on mass transport using three-dimensional modeling.

Comparison between Ocean and Murakami heat flux models
Delft3D has five different heat flux models. In this study, we analyzed two of them: Murakami and Ocean heat flux models, since both of them consider the absorption of solar radiation as a function of depth. Both models calculate the effective back radiation (Q eb ), which is the difference between long wave radiation emitted from the water and the atmospheric radiation, the evaporative heat flux (Q ev ), the convective heat flux (Q co ) and net solar radiation (Q sn ). The total heat flux (Q tot ) is: Q tot =Q sn -Q eb -Q ev -Q co . Figure 5 shows the total heat flux (Q tot ) at the air-water interface calculated for Vossoroca reservoir. When the reservoir loses heat, this loss is higher when applying the Ocean model. The mean energy flux was 102.8 Wm -2 and 10.3 Wm -2 for Murakami and Ocean model, respectively. We investigated the terms in the total heat flux equation, and the highest differences were observed in the evaporation and convection terms (Table 2). In the evaporation term, the mean heat flux was 3.32 Wm -2 and 72.3 Wm -2 for Murakami and Ocean models, respectively. In the convective term, the mean heat flux was 1.25 Wm -2 for Murakami model, and 29 Wm -2 for Ocean model. Comparing the equations for the two models and a after running sensitivity analysis, the cause of this difference was identified being a parameter in the evaporation equation. In the wind speed function: f(U 10 ) = cU 10 , c=1.2 10 -9 for Murakami model, while in Ocean model, this parameter is 0.0015. The term related to convection had the same shape and behavior as the evaporation term, because convective heat flux is calculated as a function of evaporation.
However, even changing the coefficient in the Murakami model, the results were not identical compared to the results calculated using Ocean model. Polli and Bleninger (2016) used measured temperatures of a closeby reservoir (Capivari reservoir also in Paraná State, Brazil), and estimated the heat fluxes. Comparing the results, the Ocean heat flux model agrees well with the fluxes estimated in this area, which is very similar (in terms of topography, vegetation, climate) to the area of Vossoroca reservoir. In addition, the Murakami model application for Vossoroca reservoir showed higher increases in water temperature, when compared to measurements and the Ocean model. The simulated temperature in a 14-layer model, using Murakami model, reached a mean absolute error (MAE) of 1.33±0.93 °C. Using the Ocean model, the error was 0.64±0.54 °C. Therefore, we decided to use the Ocean heat flux model for further studies to investigate thermal stratification and additional analysis in Vossoroca reservoir.

One and three-dimensional heat transport modeling
The modeling of Vossoroca reservoir was performed for a period of approximately three months. Table 3 provides general information of the Vossoroca reservoir model set-ups for the 1D and 3D models. MTCR-1 used 85 uniform fixed layers (0.2 m), a 2 min time step and its processing time is around 1 min. Delft3D, on the other hand, was set-up with 20 uniform fixed   layers (Z-model), a 3s time step (due Courant criterion), and a horizontal grid of 15 m x 15 m. The time to run Delft3D in these conditions, in parallel (with 4 processors), was around 10 days. Both models were calibrated considering measured temperature profiles at the floating platform (see Figure 1 and Figure 2). Figure 6 shows results for temperature profiles at the point of the floating platform. Both models simulate the initial period of stratification and mean temperatures. MTCR-1 (Figure 6b) tends to increase temperature more in deeper regions than Delft3D when stratification is present, but both models follow reasonably well the measured temperature profiles. The quantitative comparison showed that Delft3D (Figure 6c) represents measured data (Figure 6a) more accurately.
The quantitative comparison was done calculating the mean absolute error (MAE), and standard deviation, and is shown in Table 4. The difference between the models reflect their underlying assumptions. The values for MTCR-1 are however similar to those of other 1D models. We also estimated the errors for the temperatures of each layer (Table 5). For the 1D model, the error varies from 0.87 °C to 1.60 °C. For the 3D model, from 0.36 °C to 0.72 °C, respectively. For both models, the standard deviation is highest at the surface layer.
In that regard we investigated the vertical temperature transport, which is based on turbulent transport, parametrized by eddy diffusivity. The parameter was calculated for both models -the 1D model uses a Richardson formulation, varying with depth and being influenced by wind and the 3D model, which uses a standard k-ε model. Figure 7 shows the variation of vertical profiles of eddy diffusivity at the point of the floating platform. The overall magnitude of eddy diffusivity for the 1D model (Figure 7a) is lower than the 3D model ( Figure 7b).
On day 07/16/2012, the reservoir is fully mixed, and high eddy diffusivity coefficients were calculated by the 3D model (with a maximum of 110x10 -5 m 2 s -1 ). On the same day, the 1D model had a maximum of 7.5 10 -5 m 2 s -1 (the mean wind speed was 1.0 ms -1 ). On the other hand, on day 07/30/2012, the maximum eddy diffusivity for both models was similar -9.0 10 -4 m 2 s -1 for the 1D model, and 9.1 10 -4 m 2 s -1 for the 3D model (mean wind speed of 2.9 ms -1 ). In addition, during the stratified period, high diffusivities were calculated close to the surface by both models, however, with lower values in the 1D model. Below the thermocline, both models simulated low eddy diffusivity coefficients. The k-ε model showed a better representation of the system than the formulation using Richardson number of the 1D model. This explains also the performance difference of the 3D model in relation to the 1D model.
For an overall, integrated analysis of stratification, vertical temperature gradients between surface and bottom were calculated for the measurements and compared to model results. Figure 8 shows time series of those gradients. Both models follow the expected behavior of the measured data but, most of the time the MTCR-1 model underestimates the observed temperature gradients. We estimated stratification periods of certain temperature gradients (see Table 6). Measurements show that a temperature gradient higher than 1 °C occurred over 79 days (out of 94 days      Table 6 also shows the periods when gradients are larger than 2 °C and 3 °C, showing the same underestimations for the 1D model. The Wedderburn number (W) and Lake number (L N ) were also calculated for Vossoroca reservoir, using temperature measurements, meteorological data and modeling results. These physical indices provide information about the mixing potential (L N <1 and W<1) or stratification stability (L N >1 and W>1) of the water column. Their critical values are L N =1 and W=1. These indices were calculated using the software Lake Analyzer (READ et al., 2011) and are presented in Figure 9 and Table 6.
For the measurements (Figure 9a), in general, Wedderburn number is higher than 1, with exception of the beginning of the period, where high wind speeds and low differences in temperature in the water column are observed. In August, the surface layer starts to warm, and the increase in temperature gradients combined with low wind speed results in higher Wedderburn numbers, corresponding to strong stratification. Lake number shows the same behavior -in the beginning of the period presented, with low numbers, indicating weak stratification and allowing mixing in the water column -which is suppressed by the strong stratification after August. Table 6 shows that in 91 days Wedderburn and Lake number are higher than 1, indicating stability of the water column. For the one and three-dimensional model results (Figure 9b and Figure 9c), most of time, the results are higher than the critical number (W=1 and L N =1), indicating stability of the water column and stratification. The number of days that Wedderburn and Lake number were higher than the critical value was close for the two simulations (see Table 6) -meaning that even the one-dimensional model did simulate temperature profiles less accurately than the three-dimensional model, physical indices show the same behavior.
An additional advantage of the 3D model over the 1D model is that it provides spatial distributions of the simulated quantities. Figure 10 shows the spatial variation of regions with high or low occurrence of periods with temperature gradients being larger than a pre-defined criterion. A temperature difference between surface and bed layer of ΔT>1 °C occurs over approximately 80% of the simulated time in most of the reservoir area -the exception is close to São Joãozinho and Vossoroca rivers that are shallow and in the side arms of the reservoir, where this gradient occurs with less frequency. A ΔT>2 °C occurs over around 60% of the time in an area close to the one where ΔT>1 °C occurs. A ΔT>3 °C varies from 0% to 40% of the time near to São Joãozinho river, and the frequency increases moving to deeper regions of the reservoir. Gradients of ΔT>4 °C occur with less frequency close to the rivers, since these regions are shallower, as well in the side arms. In the deepest area of the reservoir, ΔT>4 °C occurs over around 50% of the time. A ΔT>5 °C is calculated as being less than 10% of the time in shallow areas -i.e. close to the rivers, but also in the side arms. A ΔT>6 °C occurs with lower frequency in all the reservoir area -in the deepest area, it is less than 20% of the time. Those results clearly show that the assumption of the 1D model (no horizontal temperature gradients) holds only for the deeper reservoir regions, and additionally varies over time.
Thus, in addition to the differences between the formulations of turbulent vertical mixing, the model differences are also related to the representation of horizontal temperature gradients. Those are partly included in the 1D model, but only in an integrated manner, by considering temperature differences and heat flux from colder or warmer river inflows. But the 1D model does not consider effects of differential cooling or heating in side arms, or alike, which are represented in the 3D model.
As 1D models are widely used for predicting stratification periods, but also for water quality analysis, we investigated the role of those horizontal temperature gradients on the associated substance mass transport as shown in the next chapter.

Tracer transport modeling
In this section, we analyze and quantify the importance of temperature gradients on mass transport from a river flowing into the reservoir and due to differential cooling or warming in a side   arm -considering that vertical and horizontal temperature variations may interfere in the mass transport, and using three-dimensional modeling. A conservative tracer was continuously released in São João river and in a side arm, as illustrated in Figure 11. For the river tracer modeling, a longitudinal section from the river until the deepest point of the reservoir was defined for visualizing the results. In addition, 8 cross-sections where chosen for the quantification of the advective and dispersive transport (Figure 11). For the side arm tracer modeling, a longitudinal section from the shallow to the deeper area was defined, as well as 5 cross-sections for the quantification of the advective and dispersive transport (Figure 11), respectively. For more information about the quantification of advective and dispersive transport in Delft3D, see Deltares (2014).

River tracer transport modeling
In this section, we analyze and quantify the path and distribution of a tracer continuously released in São João river (see Figure 1), considering water density effects. The same configuration and model set-up as shown in the previous sections was used (see Table 3), just including the tracer in the three-dimensional modeling. Figure 12 shows the results visualized along the longitudinal section for temperature and normalized concentration (C/C max ) for four different time steps of the simulation. As it is shown, depending on the river and water column temperatures, the path the river water flows may vary and affect reservoir dynamics in different ways. On day 07/10/2012, the river water flowing to the reservoir had a temperature of 15 °C, and the resulting tracer concentration indicates an intrusion flow occurring. At the second time step presented, day 07/18/2012, river temperature was 12 °C, a temperature which is also verified near to the bottom of the reservoir. Therefore, an underflow is occurring, as identified by the tracer distribution and illustrated by the line of 30% of maximum tracer concentration close to the bottom. On day 08/05/2012, temperature of river water was 17 °C, which is higher than surface temperature (mean temperature of 16.3 °C). Thus, an overflow occurs, as shown by the tracer distribution over the water column. During the last time shown in Figure 12, day 08/28/2012, an

11/19
intrusion is verified again. The results also show that the tracer reaches the deepest area of the reservoir after approximately 18 days travel time.
In order to analyze the effects of temperature on mass transport, we performed the same simulation as presented, but without including heat transport in the modeling. Figure 13 shows the results of normalized tracer concentration (C/C max ) over the water column. In comparison to the simulations with heat transport results show that there is no defined path of the tracer -the tracer spreads over the entire water column. The travel time necessary for the tracer reaching the deepest area of the reservoir is approximately 1 month, compared to 18 days, when including temperature in the modeling.
For Vossoroca reservoir the modeling thus showed that temperature modeling is important, since the path the water flows is a function of water column temperature. Inflows may cause water quality variations and they are relevant for the management of these system. This has also been reported in monitoring studies.
Inflows to the reservoir may act as a source of particulate material that may cause water quality variations (adding materials to the more productive zones) (MARTIN; MCCUTCHEON, 1999), increased turbulence or introducing oxygen in deeper layers. Li et al. (2015) showed that storm-rainfall events (which introduce high inflows to the reservoir) reduced stability as a function of inflow temperature. Huang et al. (2014) analyzed the path of turbid density currents in storm-rainfall events, to identify an optimal intake location for a drinking water treatment plant. These studies show how inflows may affect thermal structure and water quality in lakes and reservoirs.
To investigate the hydrodynamic effects of temperature gradients on mass transport in more detail, we analyzed modeled velocity profiles (Figure 14) at three points R1, R2 and R3 (see the points in Figure 12) at the same time steps presented previously. R1 corresponds to the location of cross-section 2 (CS 2), R2 corresponds to the location of cross-section 5 (CS 5) and R3, cross-section 7 (CS 7). When temperature is not considered in Figure 13. Normalized tracer concentration in a longitudinal section of Vossoroca reservoir (in this simulation, temperature is not modelled).
the modeling, the velocity profiles at the points R1, R2 and R3 are more uniform and no peaks were identified (Figure 14). This is different when including temperature modeling where some velocity peaks may be identified and less uniform profiles can be observed. At 07/10/2012, the velocity peak in R1 corresponds to the location of the intrusion layer of the tracer. R2 and R3 also present higher velocity near this layer. At day 07/18/2012, the underflow induces velocities near the bottom layer -R3 has bottom velocities between 0.03 and 0.04 ms -1 . At day 08/05/2012, the tracer moves in the surface layer -where R2 and R3 are also influenced by wind direction (North-east).
At 08/28/2012, there is a peak in all three profiles R1, R2 and R3 -which also corresponds to the intrusion layer, the velocity at R3 is around 0.04 ms -1 . Figure 12 and Figure 14 show that temperature has a direct influence on the dynamics of the system -and associated velocities may increase and substances may be transported according differently and faster due to the related density effects.
To quantify the overall effects of temperature on mass transport, we calculated the advective and dispersive transport of the tracer for the model with and without heat transport. The difference between the two calculations was considered as the effect of temperature on mass transport. Table 7 shows the difference in percentage on advective and dispersive transport in the simulations. A positive sign means that that mass transport Figure 14. Velocity profiles in points R1, R2 and R3 (see Figure 12) in the river tracer modeling. *Positive sign means that mass transport through the cross-section is higher in the modeling with temperature; Negative sign means that mass transport through the cross-section is higher in the modeling without temperature through the cross-section is higher for the model with temperature effects. In general, temperature effects cause increased advective transport in all cross-sections, with exception of CS 1, in which there is a decrease. The increase in cross-section CS 8, for example, was approximately 40%. For the dispersive transport, on the other hand, when temperature is included, there is a decrease -that may reach 421% difference in CS 8. CS 1 is an exception, there was an increase in dispersive transport. CS 5 and CS 7 presented no difference in dispersive transport. The differences in transport were higher in dispersive transport, but in magnitude, advective transport is more important.
The tracer modeling showed that temperature has an important role in mass transport in lakes and reservoirs -both advective and dispersive transport are affected by temperature. Advective transport increased when temperature is considered in the modeling. Dispersive transport, on the other hand, tend to decrease when temperature is considered (but it is lower in magnitude in relation to advective transport).

Side arm tracer transport modeling
In this analysis, a small flow was included (less than 1% of the mean flow) entering the reservoir in the side arm (see Figure 11), in a way that the flow velocity was not influenced by it, and a concentration of a conservative tracer was defined and continuously released. The main goal of this simulation was to identify if temperature may be responsible or enhance mass transport from or into side arms in dentritic reservoirs, such as Vossoroca. Figure 15 and Figure 16 show water temperature and normalized tracer concentration (C/C max ) every two hours at 09/23/2012 between 00h and 14h. At this day, the side arm was stratified with temperatures near to 22 °C in the surface layer, and 14 °C at the bottom. At 00h, the tracer is confined near the shallow area of the side arm, but velocity vectors plotted with temperature show that there is an increase in velocities close to the layer of -6 m, induced by temperature differences between shallow and deeper areas. Also, in the surface layer, wind speed acts in North-east direction, but velocity vectors show the water at the surface moving in East direction. At 02h, the tracer is confined in the shallow region and velocities are still higher in the layer of -6 m. Between 04h and 08h, the tracer spreads and surface temperature decreases -the graph also shows the location of where 30% of the maximum concentration is found in the water column.
At 10h, an intrusion of the tracer is identified in the layer of -6 m -velocity vectors also show higher velocities in this region -30% of the maximum concentration is found around 250 m from the point where the tracer is released. At 12h, no intrusion is identified, the tracer is confined in the shallow region -temperature in the surface layer increased, therefore, temperature did not induce further tracer movement. At 14h, higher velocities are identified near the surface layer, changing the direction of the water movement -the tracer is confined in the shallow region.
In order to analyze the effects of temperature on mass transport, we performed the same simulation as presented, but not including heat transport in the modeling. Figure 17 shows tracer concentration on 09/23/2012 between 08h and 12h. In this simulation, the tracer is confined in the shallow region and no intrusion is verified. At 10h, 30% of the maximum concentration is found around 100 m from the point where the tracer is released.
To investigate the effects of temperature on mass transport, we compared velocity profiles in two points (P1 and P2) in the side arm shown in Figure 15. Figure 18 shows the velocity profiles at points P1 and P2 in the modeling with and without temperature. In P1, the deeper point, when temperature is not simulated, the velocity profile is more uniform. In the modeling with temperature, on the other hand, the profile is not uniform, and changes in the direction where the tracer flow was identified. In P2, the velocity profile in the modeling without temperature is also more uniform -no flow is induced. In the modeling with temperature, on the other hand, some peaks may be identified in the modeling and they agree with the layer of intrusion (around -6 m), identified in Figure 15 and Figure 16. The higher velocities are around 0.02 ms -1 . Therefore, temperature may induce mass transport in side arms of dentritic reservoirs and also be a substance source for deeper waters. To quantify this effect of temperature on mass transport, we calculate the advective and dispersive transport of tracer for the modeling considering heat transport and the modeling without heat transport. The difference between the two calculations was considered the effect of temperature over mass transport. Table 8 shows the difference in percentage on advective and dispersive transport in the simulations. A positive sign means that that mass transport through the cross-section is higher in the modeling with temperature. In general, when temperature is considered in the modeling, advective transport increases in all cross-sections, with a maximum increase of 9.1%  in cross-section 1. For dispersive transport, on the other hand, no defined pattern may be defined -in some cross-sections there is an increase in transport when temperature is considered in the simulation (cross-sections 2 and 3). On the other hand, there is a decrease in dispersive transport when temperature is considered (sections 1 and 5). Advective transport showed to be more important than dispersive transport in the modeling.
These findings corroborate with monitoring results shown by Martin and McCutcheon (1999), where the pattern of flows and mixing in a lake/reservoir are affected by bathymetry, thermal structure, inflows, outflows and wind mixing. Monismith et al. (1990) measured temperature profiles in three different regions of the Wellington reservoir (in shallow and deeper regions) and the authors found out there is a non-one-dimensional response of the reservoir to heating and cooling, due to sloping bathymetry generating gravity currents (WELLS; SHERMAN, 2001). This response is due differential cooling and/or heating. The former is, for a given heat flux out of the water surface, temperature of the water column decreases more rapidly in the shallower regions where the water column has less thermal mass than the deeper regions (WELLS; SHERMAN, 2001). This produces a horizontal temperature gradient (i.e., water density) which is driven to the bottom as a cold gravity current that transports denser water and even contaminants, from the shallow to the deeper regions (WELLS; SHERMAN, 2001). The differential heating, on the other hand, occurs when solar radiation heats shallow waters faster than open waters (HORSCH; STEFAN, 1988). The effects of differential cooling/heating depend on the weather conditions and may be switched on or off within minutes or hours (IMBODEN; WÜEST, 1995). James and Barko (1991) measured temperature transects in Eau Galle reservoir (USA) and used a tracer (Rhodamine WT red fluorescent dye) to analyze horizontal transport due differential cooling. With the results of the hourly volumetric flow rates estimated from the tracer (between littoral and pellagic regions), they calculated the total phosphorus (TP) exchange rates, assuming that TP moves with the tracer. The measurements by James and Barko (1991) show that the higher the horizontal temperature gradients, higher the hourly volumetric flow rates.
The results and related literature studies thus showed clearly the importance of modeling vertical variations when working with reservoir modeling, and the need to link with temperature modeling. Interestingly, there are still several applications of lake modeling Figure 18. Velocity profiles in points P1 and P2 (see Figure 15) in the side arm tracer modeling in Vossoroca reservoir. Table 8. Difference in percentage between advective and dispersive transport considering temperature and not considering temperature in the side arm tracer modeling.
using 2D models or insufficient vertical resolutions of 1D or 3D models. Such approaches may have significant consequences on mass flow. To highlight those differences, we also performed a 2D simulation with Delft3D (1 layer in the vertical direction), and we analyzed the path and distribution of the tracer concentration, and effects of this approach on resulting mass transport for the case of tracer released in the side arm, and compared those results with 3D modeling. Figure 19 shows the normalized tracer concentrations for P1 and P2 points (see Figure 15). For the 2D modeling, the tracer concentration is uniform over the water column and no stratification effects can be modelled. The intrusion modelled by the 3D model (increase in velocity and tracer concentration) is therefore not simulated by the 2D model. In the deeper point (P1), the highest difference in the mean tracer concentration was at 09/23/2012 14h -the concentration in the 2D modeling was 30% higher than the mean concentration in the 3D modeling. At point P2, the mean tracer concentration of the 3D model is not reproduced by the 2D model -on 09/23/2012 06h, for example, the mean tracer concentration simulated by the 3D model is 140% higher than the concentration calculated by the 2D model. Therefore, the 2D modeling does not reproduce the mean concentrations that are expected to occur in the side arm modeling. No vertical variations are possible to simulate when a model with one layer in the vertical direction (2D model) is used or when the model does not consider temperature variations, and this affects directly mass transport in lakes and reservoirs.

CONCLUSION
Presented modeling approaches showed different performances when comparing measurements with simulations. One-dimensional models allow for identification of the general behavior of the reservoir in terms of defining periods when stratification occurs. The advantages of this approach are related to the fast-computational time of the simulation, together with the low data demand, providing the possibility to simulate several scenarios, with less input and calibration data, compared to 3D models. However, the results showed that the underlying assumptions, i.e. the chosen turbulence model, and the consideration of non-existing horizontal temperature gradients do affect accuracy. In contrast, three-dimensional modeling can simulate the observed data accurately, but with the disadvantage of high processing time. Comparing model performance by using commonly used physical indices instead, showed that the one-dimensional modeling may give very similar information as three-dimensional models, in terms of stability or mixing of the water column. Therefore, one-dimensional models are recommended to identify if the reservoir stratifies and for how long, as a first approach, and complemented with three-dimensional modeling when horizontal substance transport is of interest.
The analysis of temperature variations on mass transport showed that temperature induced density effects tend to increase the magnitude of velocity driving significant substance transport in the water column, especially in horizontal direction. From the simulations and quantification of transport, we identified that Figure 19. Normalized tracer concentration in points P1 and P2 (see Figure 15) in the 2D and 3D side arm tracer modeling in Vossoroca reservoir.
Polli; Bleninger 17/19 advective transport is enhanced by temperature (which is the most significant transport path calculated) -reaching up to 40% increase in a cross-section in the river tracer modeling, and 9% increase in a cross-section in the side arm tracer modeling. In the river tracer modeling, there was a decrease in dispersive transport considering temperature effects. In the case of the side arm tracer modeling, no pattern was possible to define in the dispersive transport.
Depth averaged 2D models by definition are not capable to identify stratification periods. And in addition, it was shown that they are not capable to reproduce the density driven substance transport neither, which for Vossoroca reservoir has been shown as significant.