Three-Dimensional Numerical Modeling of Large-Scale Free Surface Flows with Heat and Contaminant Transfer in Reservoirs

'e temperature distribution and pollutant distribution in large reservoirs have always been a hotspot in the field of hydraulics and environmentology, and the three-dimensional numerical modeling that can effectively simulate the interactions between the temperature fields, concentration fields, and flow fields needs to be proposed. 'e double-diffusive convection lattice Boltzmann method is coupled with a single-phase volume of fluid model for simulating heat and contaminant transfer in large-scale free surface flows.'e couplingmodel is used to simulate the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir. 'e mechanism of convection-diffusion, gravity sinking flow, and the complexity of the temperature and the pollutant redistribution process are analyzed. Good agreements between the simulated results and the reference data validate the accuracy and effectiveness of the proposed coupling model in studying free surface flows with heat and contaminant transfer. At last, the temporal and spatial variations of flow state, water temperature stratification, and pollutant transport in the up-reservoir of a pumped-storage power station are simulated and analyzed by the proposed model.'e obtained variations of the flow field agree well with the observations in the physical model test and in practical engineering. In addition, the simulated temperature field and concentration field are also consistent with the general rules, which demonstrates the feasibility of the coupling model in simulating temperature and pollutant distribution problems in realistic reservoirs and shows its good prospects in engineering application.


Introduction
A great number of large reservoirs have been built in China. By 2017, there are 732 reservoirs with a water capacity of over 100 million cubic meters [1]. e reservoirs play a crucial role in hydroelectric power, flood control, navigation, fisheries, and other comprehensive benefits. However, the water temperature stratification and pollutant accumulation exist widely in these huge reservoirs, which cause a detrimental effect on the ecological environment [2]. erefore, to study the temperature and the pollutant redistributions in reservoirs, an effective and accurate numerical model for simulating heat and contaminant transfer in large-scale free surface (LS-FS) flows needs to be proposed, which is beneficial to the harmonious hydraulic development and ecological protection [3].
At present, the general methods for studying heat and contaminant transfer in LS-FS flows are mainly two-dimensional (2D) numerical models, especially shallow water equations coupling convection-diffusion equations for simulating planar flow problems [4,5]. e 2D numerical simulations have been one of the most commonly used methods to analyze large-scale flows with heat and contaminant transfer, such as flows in reservoirs [6], lakes [7], and oceans [8]. However, the 2D method could not analyze the vertical temperature and pollutant distribution, and the plane simulation results are also not accurate enough because of the lack of vertical velocity effect, which evidently limits the applications of 2D models in engineering practice. e flows with an obvious three-dimensional (3D) convection-diffusion mechanism shall be simulated by more promising 3D methods. Typically, 3D numerical methods for simulating LS-FS flows are mainly based on gas-liquid two-phase flow models [9], interface tracking models [10][11][12], and interface capturing models [13][14][15]. Nevertheless, the flow problems in the field of hydraulic projects are usually complicated because of the large-scale and complex boundary conditions, which not only results in a dramatic amount of computing resource consumption but also makes the stable and accurate simulation challengeable [16].
Single-phase free surface lattice Boltzmann (SPFS-LB) model was originally proposed by ürey in 2003 [17]. Compared with the traditional two-phase flow models in which both gas and water are simulated, the SPFS-LB model is more powerful to simulate 3D LS-FS flows, because it only simulates the water flows and ignores the movement of air, so the consumption of computational resources is lower [17,18]. e lattice Boltzmann models for solving convection-diffusion problems can be mainly classified in the multidistribution function approach [19], the multispeed approach [20], and the combined approach [21]. Among the above approaches, the multidistribution function approach is the most mature and has been used successfully in many fields [19] but has seldom been applied to heat and contaminant transfer in LS-FS flows.
To study the 3D temperature and pollutant distribution in large reservoirs, the double-diffusive convection LB method is coupled with the SPFS-LB model. e rest of the paper is organized as follows. Firstly, the basic ideas and algorithms of methods are briefly described, and the treatments of the coupling scheme for improving stability and precision are addressed. en, the accuracy and effectiveness of the proposed method are verified by two benchmark cases: the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir. At last, as an attempt of simulating engineering practice, the process of the flow state, water temperature, and pollutant variations in a practical reservoir is simulated and analyzed.

e Single-Phase Free Surface LB Model.
Unlike traditional two-phase flow methods that set up different distribution functions (DFs) for water and air, the SPFS-LB model only simulates water flows (the denser and more viscous phase) by the single-phase LB equations and tracks the free surface by a constructed boundary treatment. e following assumptions should be made before the simulation by the SPFS-LB model [17]: (1) the effects of air flows (lower density and viscosity) on water flows can be neglected; (2) the air can reach the equilibrium state immediately after the state of water is changed; (3) the simulated water is separated by a closed layer of the interface cells (the free surface) from air. Absolutely, the large-scale water flows in reservoirs are in line with the above assumptions. ree types of computational cells are introduced to capture the interface, including the empty cells, interface cells, and filled cells, as shown in Figure 1. e empty cells do not contain any water, so there is no need to define physical quantity. e filled cells are full of simulated water, in which the definition of physical quantity and the evolution processes are the same as those in the traditional single-phase LB models. e interface cells are tracked by computing the water volume fraction ε, which can be defined as ε � m/ρ, with m and ρ being the mass and density of cells, respectively. And the values of empty cells, interface cells, and filled cells are 0, 0-1, and 1, respectively.

Computation of Mass Flux.
e movement of the interface in the SPFS-LB model is realized by the transformation of cell types. And the cell type is judged by the volume fraction ε which can be updated by calculating the inflow and outflow of cell mass. By conducting discretization in time and on uniform lattices and applying the divergence theorem to the convective term, the continuity equation in discrete form can be derived as m ( where Δm i denotes the mass variation on ith direction along with the particle velocity [22]. In the LB method, the mass flux can be simply computed through the two antiparallel particle distribution functions f i and f inv(i) (see equations (1) and (2)), where e inv(i) � −e i [23]. 1, x + e i ∈ empty cell, interface cell, filled cell.
erefore, the mass m of the interface cell at lattice node x for the time step t + Δt can be updated the following equation, which indicates that the updated m is obtained by summing the mass flux in all velocity directions.
where n denotes the number of discrete velocities. It can be proved that the mass updated by equation (3) is conserved [22].

Reconstruction of Interface Cells.
In the SPFS-LB model, the propagation and collision of LB equations can be normally performed in filled cells for that they are not adjacent to empty cells. However, the interface cells are always surrounded by empty cells, where physical quantity and DFs are not defined (shown in Figure 1). According to assumption 2, the missing DFs of the interface cells can be reconstructed based on the macroscopic boundary conditions and the force balance for opposite lattice directions [17]. e following equation shows the reconstructed DFs of interface cells at x with surrounding empty cells at x + e i Δt, where u and ρ G denote the velocity and the pressure of interface cells at position x, respectively.
It should be pointed out that not only the missing DFs but also all the DFs whose discrete velocity e i satisfies e i n < 0 need to be reconstructed through equation (4) [17]. e vector n refers to the surface normal direction, which is obtained by the second-order central difference approximation n � s i e i ε(x + e i ).

Mass Allocation.
Interface cells need special treatments after the updating of their mass m and volume fraction ε. e interface cells whose updated ε is greater than 1 shall be converted to filled cells, and the excessive mass needs to be allocated to the adjacent interface cells and empty cells. After getting mass, the surrounding empty cells shall be turned into interface cells, and their macroscopic variables can be obtained by averaged velocity u and density ρ of the surrounding nonempty cells.
e DFs of new emerging interface cells can be initialized by the equilibrium distribution functions of the LB method [18]. Similarly, the interface cells whose ε is less than 0 shall be converted to empty cells and the surrounding filled cells turn to interface cells. e negative m needs to be compensated to 0 by the neighboring nonempty cells. e proportion of the excess mass distribution and negative mass compensation can be obtained by the dot product between the interface normal direction and the distribution direction. e equations can be expressed as [17] where m ex and υ refer to the excess mass and proportion, respectively. e interface moves as the completion of mass allocation and cell type conversion.

Double-Diffusive Convection Lattice Boltzmann Method.
From a large-scale point of view, the water flows in reservoirs are normally considered incompressible, whose viscosity, density, thermal diffusivity, and contaminant diffusivity can also be assumed to be constant except for the buoyant. erefore, the governing equations for simulating heat and contaminant transfer in reservoirs can be expressed with the Boussinesq assumption [21] as where F, ], α, and D denote body force term, viscosity, thermal diffusivity, and contaminant diffusivity, respectively. In [24], the double-diffusive convection lattice Boltzmann (DDC-LB) method is used to solve equation (6). e DDC-LB method simulates velocity, temperature, and concentration fields by three sets of respective DFs, and the model couples the three fields by a body forcing in the velocity field. Here, the method is extended to 3D form on the basis of the 2D DDC-LB method in [24], which simulates velocity field by D3Q19 model and uses D3Q6 model for temperature fields and concentration fields simulation. e LB equation of the 3D model can be described as i (x, t) , i � 1, 2, 3, 4, 5, 6, i (x, t) , i � 1, 2, 3, 4, 5, 6. Advances in Civil Engineering f i , g i , and h i refer to the DFs for the velocity field, temperature field, and concentration field. Variable τ denotes the relaxation factor, which can be determined by the viscosity and diffusivity, and Prandtl number (Pr) can be obtained by τ (see the following equations): e equilibrium distribution functions f eq i , g eq i , and h eq i in equations (10)- (12) can be expressed as To derivate governing equation (6) by Chapman-Enskog expansion, coefficients σ, λ, and c need to satisfy λ + 2c � σ and λ + 2c � 1/3, which are chosen as 1/3, 1/6, and 1/12 in this work, respectively. s i (u) is a function of the discrete velocity u: ω i is the weight coefficient, ω 0 � 1/3, ω 1∼6 � 1/18, and ω 7∼18 � 1/36 in the D3Q19 model. e macroscopic quantities u, P, T, and C are the macroscopic velocity, pressure, temperature, and concentration; they can be obtained by the first or the secondary moment of the corresponding DFs as Simulating 3D temperature fields and concentration fields in large-scale flows by the D3Q6 model can reduce calculation amount and save computer memory usage. However, when precisely analyzing the small-scale flow problems with complex boundary conditions, using the D3Q19 model instead of the D3Q6 model can retain better accuracy and stability. e number of computational grids for simulating the model reservoir (in Section 3.2) is over 10 million (2440 × 46 × 92), and there are over 37 million (1024 × 1024 × 36) for simulating the practical reservoir (in Section 4). e computer memory usage and computation time are estimated to be increased by 45% and 23%, respectively, if the temperature fields and concentration fields are simulated by the D3Q19 model instead of by the D3Q6 model. Due to the constraint of computer performance, the D3Q6 model is used for the present work. Good agreements between the simulated results and the reference data in Section 3 show that the results simulated by the D3Q6 model are accurate enough.

Key Treatments of the Coupling Scheme.
e SPFS-LB model is firstly coupled to the DDC-LB method in the present work, and several key treatments need to be adopted to make the coupling scheme compatible and to obtain reasonable simulated results.

Reconstruction and Initialization of Temperature Fields
and Concentration Fields in Interface Cells. As described in Section 2.1.2, the undefined DFs of empty cells would enter into interface cells, so the DFs of temperature fields and concentration fields in interface cells need to be reconstructed after the LB propagation. According to assumption 2 (introduced in Section 2.1) and referring to the reconstruction of flow fields, the missing DFs in the temperature field and concentration field can be calculated by the following equations: 4 Advances in Civil Engineering Some empty cells would convert to interface cells after the update of cell types, and the temperature T, concentration C, DFs g i , and DFs h i should be initialized. is coupling model neglects the exchange (absorption or release) of heat and contaminant between water and air. erefore, the macroscopic T and C of new emerging interface cells originate from the surrounding excessive cells, and they can be obtained by the weighted average of the variables of surrounding new filled cells. e weight can be computed according to the proportion of the assigned excess mass. Referring to the initialization of DFs f i , DFs g i , and DFs h i in new converted interface cells can be initialized by equations (10)-(12).

2.3.2.
e Force Term for the Buoyancy Flows. e flow velocity dominates the transfer of heat and contaminant, while the buoyancy induced by the nonuniform distribution of temperature and concentration has an impact on the water flows. To simulate the buoyancy flows and realize the two-way coupling between the three fields, an additional body forcing term Fi needs to be attached to equation (7). e present model adopts the approach proposed by Cheng and Li [25] on account of its stability and accuracy for treating nonuniform, unsteady source terms. e forcing term Fi can be calculated by the following equation according to [25]: where in which A and B denote the source term in the continuity equation and momentum equation, respectively. Considering equation (6), the present model sets

Large Eddy Simulation for Turbulence.
Turbulence exists in practical large-scale water flows with heat and contaminant transfer. e subgrid-scale model based on large eddy simulation is introduced to the coupling scheme. e physical variables are separated into large-scale variables and small-scale variables by a certain filter function G in the subgrid-scale model [26,27]. Take φ � φ + φ′, for example. e large-scale variable φ can be expressed as follows. φ′ presents the small-scale pulsation, which is modeled in a simplified way. e following equation should be obtained through filtering equation (7) and assuming f e relaxation time τ e is not a constant. It varies with space and time, which can be calculated by ] e through equation (8). e variable ] e can be understood as the total viscosity ] e � ] 1 + ] t , in which ] 1 and ] t denote the molecular viscosity and eddy viscosity (or turbulence viscosity), respectively. ] 1 is determined by the physical property of water, and ] t can be computed by ] t � (CΔ) 2 |S|, where C, Δ, and S refer to Smagorinsky constant, filter width, and strain-rate tensor, respectively. According to [27], the strain-rate tensor S can be easily computed by using the second-order moments of the filtered particle distributions.

2.3.4.
e Coupling Procedure of the Coupling Scheme. is section makes a summary of the calculation procedure of the coupling scheme proposed in the paper. Before the iteration, the geometry needs to be initialized, and the reasonable boundary condition should be adopted in accordance with the simulated flow problems. en given all the variables at time step n, the procedure for marching to step n + 1 is as follows: (i) Perform LB propagation with the boundary conditions (ii) Compute the mass flux, and update the mass of cells as described in Section 2.1.1 (iii) Calculate the volume fraction of nonempty cells (iv) Reconstruct the DFs of the flow field, temperature field, and concentration field by equations (4) and (15)

Model Verification
To validate the feasibility and accuracy of the proposed coupling scheme, the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir are simulated, and the detailed comparison and analysis are made between the simulated results and the reference data. If not otherwise stated, the acceleration of gravity, water density, and kinematic viscosity are set as g � −9.8 m/s 2 , ρ � 10 3 kg/m 3 , and ] � 10 −6 m 2 /s, respectively.

Double-Diffusive Natural Convection in a Cubic Cavity.
In the heat and contaminant transfer flows, the doublediffusive natural convection in a cubic cavity reflects the essential problem of two-way coupling between temperature field/concentration field and flow field. At present, there have been many classic physical experiments or numerical research results [28,29]. In a cubic cavity of L × L × L, the front (wall 5) and back (wall 6) are smooth walls, and the rest (wall 1∼4) are nonslip walls (see Figure 2). e temperature and concentration at the left and right boundaries are kept constant at T 0 /C 0 and T 1 /C 1 , respectively (T 0 < T 1 ; C 0 < C 1 ); the top and bottom are adiabatic for heat transfer and impermeable for contaminant transfer; that is, the derivatives of temperature fields and concentration fields are always 0 along the normal direction of the boundary. Initially, the cubic cavity is filled with homogeneous fluid with velocity u � 0, temperature T � (T 0 + T 1 )/2, concentration C � (C 0 + C 1 )/2, and Pr � 0.71. e flow characteristics of double-diffusive natural convection can be estimated by the following three dimensionless numbers: temperature Rayleigh number Ra T � (gβ T (T 0 − T 1 )L 3 /υα), concentration Rayleigh number Ra C � (gβ C (C 0 − C 1 )L 3 /Dα), and Lewis number Le � (α/D). e initial values of temperature T 0 /T 1 and concentration C 0 /C 1 can be assigned to arbitrary numbers that satisfy the formulas above, such as T 0 /C 0 � −1 and T 1 /C 1 � 1 in the present simulation.
In order to quantify the simulated results, several dimensionless numbers along the right hot wall are calculated below, which are local Nusselt number, average Nusselt number, local Sherwood number, and average Sherwood number, respectively, as shown in the following equations: Firstly, the Rayleigh number and Lewis number are set as Ra T �10 4 , Ra C � 10 4 , and Le � 1, respectively. e fluid in the cavity is driven by both temperature field and concentration field. e added driving flows are generated if the driving directions of temperature and concentration are the same.
Otherwise, the opposed driving flows are generated. Both flows are simulated in this section. e computational grid is set as 128 × 128 × 128. It should be noted that the simulated distribution of temperature field and concentration field is exactly the same, for the same convection-diffusion equation is applied to simulate heat transfer and contaminant transfer, and the same calculation parameters are set in the present simulation in order to compare with the research results in [29]. Figures 3 and 4, respectively, show the streamline and isothermal surfaces (or isoconcentration surfaces) after the convergence of simulation. For the added driving flows, the isothermal surfaces move farther on the upper and lower wall compared with the single-diffusive convection under the same Rayleigh number, which shows that the heat and contaminant spread more fully in the cubic cavity. e reason is that the flow field is driven by the buoyancy from temperature difference and concentration difference simultaneously, and the convection that leads to the spread of heat and contaminant is more significant. As shown in Figures 3(a) and 4, driven by the downward buoyancy around the left wall and the upward buoyancy around the right wall, the thin boundary layers are formed on the two walls. e flow forms a stable cycle, while a stagnation area appears in the center of the cubic cavity.
For the opposed driving flows, the buoyancy caused by temperature difference and concentration difference is exactly equal in magnitude and opposite in direction, and their driving effect on the flow field also cancels each other out. e fluid in the cavity does not flow at all (see Figure 3(b)). In this case, the propagation of heat and contaminant is completely dependent on diffusion, which can be proved by the distribution of equidistant vertical isothermal surfaces/

Non-slip Wall
Wall (3) Smooth Wall Wall (6) Wall (2) Wall (5) Wall (1) Wall (  Advances in Civil Engineering isoconcentration surfaces in Figure 3(b). e above results are consistent with the research in [28,29]. e local Nusselt number of the added driving flows along the right wall (hot wall) is calculated and plotted in Figure 5. In the region of Z/L � 0∼0.13, the isothermal surfaces are vertically distributed, and the temperature gradient along the Y direction is basically unchanged (see Figure 3(a)). e Nu rises slightly from 4.18 to 4.43 with the increase of boundary layer velocity. After reaching Z/ L � 0.13∼0.96 regions, the velocity increases sharply and the transverse convection enhances, leading to a significant decrease of the temperature negative gradient −zT/zy along the Y direction. erefore, Nu plummets from 4.43 to 0.64. Subsequently, influenced by the upper wall, the velocity and temperature gradient remain relatively stable, and Nu also stabilizes at about 0.64. It can be seen from Figure 5 that the simulated results are very close to the data in [29], except for slight differences on the upper and lower corner of the right wall (the maximum relative error is 0.72%). e relative error is inevitable, and it is enlarged especially at the corner, for different numerical methods and boundary treatments are adopted in this paper and the reference.
Additionally, the added driving double-diffusive convection with the temperature Rayleigh Ra T �10 7 and concentration Rayleigh number Ra C � 10 6 , 5 × 10 6 , 1.5 × 10 7 , and 5 × 10 7 is simulated. e average Nusselt number and average Sherwood number are calculated along the hot wall. According to Table 1, the dimensionless numbers from the present simulations agree well with the results in [29], and the errors are all within 1% (the literature results only kept a decimal place, so it is difficult to measure the errors more accurately). ese analyses demonstrate the accuracy of the

Temperature Distribution in a Model
Reservoir. In this section, the variation process of temperature field and flow field in a model reservoir is simulated by the proposed coupling model, so as to demonstrate its effectiveness in simulating practical physical problems. e simulation is based on the data of gravity undercurrent experiments, tested by the US Army Corps of Engineers, Waterways Experiment Station [30,31]. e geometry is shown in Figure 6. e inlet is a 0.15 m high hole at the bottom of the left-hand wall, while the outlet is a banded hole that is 0.04 m high and 0.15 m away from the bottom of the right-hand wall. Initially, the reservoir water is stationary with a uniform temperature of 21.44°C. After the simulation starts, the 16.67°C cold water is discharged at a flow rate of 0.00063 m 3 /s, while the water flows out from the outlet at the same flow rate.
Due to the lateral symmetry of the simulation, only a half is selected for calculation and results show (see Figure 6). Correspondingly, the computing grid is 2440 × 46 × 92, and the symmetric plane is treated as the symmetric boundary condition. e inlet and outlet are set as a velocity boundary condition of 0.014 m/s and 0.0173 m/s, respectively, and the rest are adiabatic and nonslipping solid wall boundaries.   (3) Horizontal thermal stratification stage: when the tongue-shaped cold current comes to the right wall of the reservoir, the thermal stratification formally formed. e backflow area at the top of the reservoir stretches across the entire reservoir and reaches the longest, as shown in Figure 8(e). As the cold inflow water forward flows out of the reservoir, its disturbance to the upper warm water decreases and the surface backflow intensity declines gradually. e major heat-transfer mechanism changes from convection to diffusion, which is also the reason why the water temperature in the reservoir ultimately presents horizontal stratification. At last (T � 50 min), the 20°C isothermal surface expands to the whole reservoir, and the 18°C and 19°C isothermal surfaces are also in parallel with the bottom. e flow condition of the whole reservoir is composed of warm backflow at the upper part and cold current at the lower part, as shown in Figures 7(f ) and 8(f ).  Based on the above analyses, the flow field and the temperature field are strongly coupled. e accurate simulation of buoyant flows is the key to study the thermal stratification of the reservoir water. Caused by the temperature difference, the inflow water sinks and the transfer of vertical heat is restricted during the whole moving process of the cold current, which strongly affects the formation, development, and stabilization of flow state and temperature distribution. It also explains the mechanism of horizontal thermal stratification in reservoirs. Figure 9 compares the velocity profile along a vertical line at 11.43 m from the reservoir inlet at t � 11 min (see Figure 6 for the position of the line), and Figure 10 shows the histories of the discharge water temperature. e maximum velocity of the cold current and upper warm water is about 0.024 m/s and −0.0065 m/s, respectively. By comparing the results simulated by the proposed coupling model with the experimental results [30,31] and the results obtained by Fluent software [32], the current layer simulated in this paper is thicker, and the velocity of backflow is smaller. Some errors exist because it cannot sensitively capture the sharp velocity gradient of the boundary layer. Besides, no boundary layer model is adopted to simulate the turbulence near the reservoir bottom. In terms of the temperature of the outflow water, the results from the present coupling model are superior to the results calculated by commercial software [32]. In summary, the above simulation shows that coupling the DDC-LB method and SPFS-LB model can simulate the large-scale thermal flows in model reservoirs effectively.

Temporal and Spatial Variations of Temperature and Concentration in a Practical Reservoir
An upper reservoir of a pumped-storage power station is about 700 m in length and width, and the maximum depth is about 25 m (shown in Figure 11). e inlet, located at the bottom of the right-hand wall, is a rectangular region within the coordinate range of X �  Figure 12 shows the streamlines in the reservoir and the velocity distribution on the slice of Z � 7 m after pumping for 5 min, 15 min, 85 min, and 100 min, among which the blue area is the water area.

Variations of Flow Fields.
Under the pumping operation modes, the water bypasses Mountain A and Mountain B ( Figure 11) and rushes to the deepest part of the reservoir basin once it flows into the reservoir. Under the influence of the bottom slope, the water flow keeps a high velocity, up to 4 m/s (Figure 12 Measured [27] RNG [29] Realizable [29] RSM [29] Present Measured [27] RNG [29] Realizable [29] RSM [29] Present

12
Advances in Civil Engineering will scour the reservoir basin and disturb the deposited sediment. As shown in Figure 12(b), the water flow at the inlet is deflected by the right bank slop of the reservoir, and it impacts Mountain A from the left side, which will cause bank erosion to Mountain A and thus affect its stability. e water flows into the open area of the reservoir basin after passing through Mountain A, forming a clockwise circulation with a diameter of over 100 m. After 15 minutes, the inflow water covers the whole basin, and the water surface calms down and gradually rises. e incoming water accounts for 65% of the reservoir capacity at the time of t � 85 min. As the flow rate at the right side of Mountain A is much larger than that at the left side, a counterclockwise circulation across the entire reservoir forms with a diameter of over 300 m, as shown in Figure 12(a). Under the disturbance of Mountain A and Mountain B, a small counterclockwise circulation (about 50 m in diameter) appears between the two mountains, and this place will be a gathering area of floating objects (such as garbage, vegetation, and others). e flow velocity decreases gradually with the rise of the water level. At t � 100 min, the flow state of recirculation zones becomes stable, and the streamlines become smooth; see Figure 12(d).
e above-obtained variation laws of the flow field agree with the observations in the physical model test and in practical engineering [33]. Figure 13 shows the isoconcentration surfaces (C � 0.4, 0.2, 0.1, and 0.05) in the reservoir after pumping for 5, 15, 85, and 100 minutes. e point pollution source is located between the two mountains. According to the analysis of flow fields in Section 4.1, the water flows from right to left at a fast speed under the action of Mountain A, and the pollutants also rapidly spread through the water before the time of t � 15 min. e spread speed of isoconcentration surface (C � 0.05) is equal to the flow velocity, and the pollutants almost do not propagate in the reverse flow direction (Figures 11(a) Advances in Civil Engineering 13 and spread to the deepest part of the reservoir basin, while their propagation distance in the clockwise direction is relatively close due to the significant inhibition effect of circulation, as shown in Figures 11(c) and 13(c). In order to control the transfer pollution, effective engineering measures need to be taken to stop the counterclockwise circulation of reservoir water. As the flow velocity in the reservoir slows down, the convection that causes the spread of pollutants gradually weakens and the diffusion effect on pollutants propagation appears. Finally, at 100 min after pumping, the pollutants almost diffuse to the entire bottom area, as shown in Figure 13(d). e above analysis of the pollutant propagation process is of certain guiding significance to the formulation of remedial measures and emergency treatment plans for pollution leakage accidents in reservoirs.  Figure 14, in which the light blue area is the water surface.

Variations of
When the cold water just enters the reservoir (t � 71 min), it sinks to the bottom due to its relatively high density. e cold water is not fully mixed with the warm water in the reservoir, and a large vertical temperature gradient is maintained at the inlet (Figure 14(a)). Five minutes later, the cold water forward reaches Mountain B, forming an undercurrent that slides into the depth of the reservoir basin at a speed of 1 m/s along the bottom slope of the reservoir. As the vertical temperature difference inhibits the transfer of vertical momentum and heat, the isothermal surfaces present a flat tongue shape, as shown in Figure 14(b). Subsequently, the advance speed of the undercurrent slows down as the bottom slope of the reservoir decreases. As shown in Figure 14(c), cold water reaches the opposite reservoir bank at t � 85 min. At the time of t � 100 min, the cold water bypasses Mountain A and reaches the deepest part of the reservoir basin. In this process, due to the slowing down of the undercurrent velocity, the dominated heat-transfer mechanism changes from convection to diffusion, and the isothermal surfaces also develop from tongue shape to pie shape, as shown in Figure 14(d). e development process and formation mechanism of vertical water temperature distribution in this reservoir are basically consistent with the analysis results in Section 3.2.
is paper summarizes the general conclusion that explains the spatial-temporal variation process of water temperature in the case of the injection of cold water into a warm reservoir. On the whole, the coupling model proposed in this paper is effective and feasible to simulate the large-scale free surface flows with the heat and contaminant transfer in reservoirs.

Conclusion
e SPFS-LB model is firstly coupled to the DDC-LB method in this paper to simulate large-scale flows with heat and contaminant transfer. To make the coupling scheme compatible and to obtain reasonable results, several key treatments described in Section 2.3 need to be adopted. e feasibility and accuracy of the proposed model are validated by two typical examples, namely, the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir. Finally, as an attempt at the field of engineering applications, the coupling scheme is applied to study the temporal and spatial variations of temperature and concentration in a realistic up-reservoir of a pumped-storage power station in detail, and some general conclusions of heat and contaminant transfer in reservoirs are drawn.
Compared with the traditional 2D numerical models, 3D temperature, pollutant, and velocity distribution can be obtained in this work. In addition, the multifield coupling mechanism of this model makes the simulation more accurate.
ese are the main advantages of the proposed coupling model. It should be noted that the proposed model ignores the exchange (absorption or release) of heat and contaminant between water and air. e nonpoint source term will be introduced in equation (7) to solve the problems in the following work. e wind and drainage also have an effect on the flow field of reservoirs, and they should be taken into consideration in the future. Besides, all the numerical simulations in the present work are performed on uniform mesh, and the practical applications require a large number of grid nodes to reach acceptable accuracy. erefore, in order to enhance the computing efficiency, other works in the future will focus on local mesh refinement to this coupling scheme and realize massive parallel computing.

Data Availability
Some or all data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
e authors declare no conflicts of interest.

Authors' Contributions
Wei Diao contributed to methodology, investigation, writing of the original draft; Hao Yuan contributed to resources and validation; Cunze Zhang contributed to funding acquisition; Liang Chen contributed to software; Xujin Zhang contributed to supervision and writing, review, and editing.