CFD model to characterize the physical processes of energy transfer in a secondary facultative lagoon

Facultative lagoons (FLs) are important technologies for wastewater (WW) treatment (Abbas et al. 2006); their importance is due to their low costs, high efficiency, simplicity of design, sustainability and easy operation and maintenance (Rao 2007; Muga & Mihelcic 2008; Olukanni & Ducoste 2011). They are very useful in developing countries, especially those located in tropical regions where there is no shortage of solar energy at any time of the year (Sah et al. 2011). The operation of an FL is affected by environmental conditions in the area where it is located. Factors such as ambient temperature, wind direction and speed, and solar radiation play an important role in the workings of this type of ecosystem (Ukpong 2013a). Therefore, the operation and design of an FL require a deep knowledge of the effects of these factors, primarily on the phenomena of energy transfer and solar radiation penetration of the body of water. These phenomena form the pillars of the ecological assemblage underlying the lagoon's decontamination processes such as nitrification, ammonification, denitrification, methanogenesis, and photosynthesis (Kayombo et al. 2000; Mayo & Abbas 2014). A limited number of studies on the energy phenomena in lagoon systems for wastewater treatment are available in the literature (Sah et al. 2012); these are related with temperature patterns of wastewater (Sweeney et al. 2005; Escalas-Cañellas et al. 2008), effect of wind on the wastewater temperature profiles (Sah et al. 2011) and light penetration in wastewater (Curtis et al. 1994; Davies-Colley et al. 2005). Studies that show the penetration of solar radiation into the body of water, the transfer of heat through the lateral walls, and the effect of the wind on these phenomena have not been reported in the literature due to their complexity and the limited experimental techniques available to characterize them. One of the techniques most recently used to study these phenomena is computational fluid dynamics (CFD) applied to lagoon systems (Karteris et al. 2005; Wu 2013). The models that use CFD have become a powerful and innovative tool to study the thermal stratification in this type of water body (Sweeney et al. 2005; Escalas-Cañellas et al. 2008). Sah et al. (2012) summarized the multiple contributions of CFD models as they applied to gaps in the treatment of WW. Since 1992, 44 CFD models have been developed for the study of FLs; of these, only 16% include some energy components (Salter et al. 2000; Sweeney et al. 2005; Sah et al. 2011). None of the models account for parameters such as the penetration of incident solar radiation or the transfer of energy through the lateral walls. However, in using these models, CFD demonstrated it has sufficient potential to explore the energy transfer processes in FLs in detail. In this context, this article contributes to the study of the physical processes of energy transfer in an FL. For this purpose, a CFD model was constructed in three dimensions (3D). Three scenarios were compared, two in which the material type in the surrounding walls was varied and another one in which the direction and magnitude of the wind on the surface were included. The model was validated with an experimental tracer study, experimental temperature data and measurements of solar radiation penetration. The results showed that the model correctly simulated the hydrodynamic of the wastewater, the temperature profiles and the penetration of incident solar radiation.

This research was performed in an experimental secondary FL in conventional configuration; it was located in the Research Station of Wastewater and Reuse in the municipality of Ginebra-Valle del Cauca, Colombia, at 3° 43′ 50″ N and 76° 16′ 20″ E, and 1,040 masl. The average temperature in this municipality is 23 °C, and the average annual rainfall is 1,280 mm, corresponding to group A (tropical) according to Köppen classification.

The inflow of wastewater into the FL was 23.76 m³·d⁻¹. The inflow was regulated using an automatic globe valve manufactured by KSB SE & Co., Frankenthal, Germany. The inlet and outlet structures are submerged tubes 0.0546 m in diameter, and they are located parallel to the longest side of the lagoon, a slope (l:z v:h) 1:1. The surface area of the ecosystem is 83.22 m², with a depth of 1.48 m, a width of 5.70 m, a length of 14.60 m, and a retention time of 3.99 days. The FL was designed for an organic load of 279 kg BOD·ha⁻¹·d⁻¹. A schematic of the experimental system is shown in Figure 1.

To validate the hydrodynamic field, a tracer study was conducted following the methodology developed by Aponte (2013). Rhodamine WT was used as a tracer using the stimulus-response technique. For this purpose, 3.8 g of the tracer was measured into a 1,000 ml volumetric flask and made up to the mark with wastewater of the FL to the tracer solution maintained the same temperature of the wastewater (Broughton & Shilton 2012). The resulting solution was added to the FL inlet in instantaneous pulse form. The tracer concentration was measured online in the effluent using fluorescence as the detection technique. For this purpose, a Turner M 8000-010 fluorometer manufactured by Turner Designs Company, Sunnyvale, California, USA, with a linear detection range between 0.4 μg·kg⁻¹ and 300 μg·kg⁻¹, was used. The wavelengths used were 550–570 nm; with the data, the distribution curves of the residence times were constructed. In the CFD model, the tracer injection was implemented comparing two methods: the discrete phase model (DPM) and chemical species model. The DPM required a higher computational expense and results were not as expected for this simulation. In contrast, the chemical species model required 80% less computational resources and the results provided were satisfactory compared with experimental tracer tests. The pulse method was used, and the concentration at the output was monitored with an ‘area weighted averaged concentration’ monitor (ANSYS 2013).

The temperature profile for the model validation was determined as suggested by Tay et al. (2012). For this purpose, three sensors were placed within the FL, and one was placed in the effluent. The three sensors were located at point P₂ (L/2) and its respective depths were 0.05 m, 0.45 m and 1.40 m. The sensors were of the reference DS1722 digital thermometer manufactured by Dallas Semiconductor, Dallas, Texas, USA, with the capacity to measure temperatures in the range from −55 °C to 120 °C with an accuracy of ±1 °C and a resolution of 0.0625 °C. Figure 2 shows the location of the sensors.

The normality of the simulated and experimental data series was validated applying the Shapiro–Wilk test (n > 50) to evaluate the normality of the data series; a P value greater than 0.05 indicates the normality. Additionally, the variances of the residence time distribution curves were compared using the Levene test (n > 50). In this case, a P value less than 0.05 indicates that there are no statistically significant differences between the variances. Accuracy tests were applied; these included the calculation of the absolute mean deviation (AMD%) and the relative mean deviation (RMD%) as measures of the error between the simulated and experimental temperature data, and the sum of the squared errors (SSE), the root mean square error (RMSE) and root mean square deviation (RMSD) as an accuracy test in the temperature data series.

A three-dimensional CFD model was developed using ANSYS Inc. Fluent^® software (V.16.1) on a Dell Precision TX 3500 workstation with an Intel Xeon X3470 processor (8 MB cache, 2.93 GHz, Turbo, HT). The 3D geometry was built using ANSYS Design Modeler^® software by ANSYS Inc. version 16.1. It was designed according to the real dimensions of the FL.

The finite volume method was used for discretization. The computational domain was divided into 161,890 hexagonal elements of 0.05 m, and ANSYS Inc. ICEM CFD™ meshing software (V.16.0) was used. The quality of the mesh elements was evaluated using the determinant and internal angle methods. The first method guarantees an element quality of greater than 0.5, and the second, an internal angle greater than 9°. This procedure was performed to favor solution accuracy and model convergence. Figure 3 shows the geometry and mesh of the CFD model.

The ‘velocity inlet’ boundary condition was used for the inlet, with an input velocity of 0.115 m·s⁻¹ corresponding to the actual fluid inlet velocity and k and epsilon parameters of 5.78 × 10⁻⁵ and 8.00 × 10⁻⁶, respectively. The suspended solid inlet concentration was 0.041% w/v. This value corresponded to the average of the experimental data measured over a period of 2 years. Two approximations were used to simulate the solids: the ‘Eulerian’ and ‘Mixture’ models. Results showed the ‘Mixture’ model decreased computing time and computational spending 13%; therefore, this model was used for the simulations; for this, the ‘pressure outlet’ boundary condition was used at the outflow. This boundary condition consists of computing the mass flow rate across the upstream zone (pressure outlet), adjusting the mass flux profile. This adjustment occurs every iteration, ensuring rigorous conservation of mass flow throughout the course of the calculation. For the outflow, k and epsilon parameters of 8.00 × 10⁻⁷ and 1.00 × 10⁻⁸ respectively, were calculated.

For the walls, the wall–temperature boundary was used to simulate the solid boundary conditions for viscous fluids and the exchange of heat between the fluid and the walls. The mixed boundary condition was used for the surface. Finally, two types of material were simulated for the lateral walls: soil and a low-density polyethylene coating (geomembrane) that covered the internal area of the lagoon.

Physical properties of the wastewater, the discrete phase (suspended solids) and the materials of the lateral walls were determined. The wastewater was assumed to be incompressible, exhibiting Newtonian behavior. The following values were assigned to the fluid: density (ρ) 1.020 kg·m⁻³ and heat capacity (Cp) 4.094 J·kg⁻¹·K⁻¹ (Stafford et al. 2014). The thermal conductivity (λ) and fluid viscosity (μ) values were 0.44 W·m⁻¹·K⁻¹ (USEPA 1995; Stafford et al. 2014) and 0.0011 kg·m⁻¹·s⁻¹ (Peterson et al. 2000) respectively. The properties of the suspended solids included a density (ρ) of 1.170 kg·m⁻³ (Huggins et al. 2004) while the particle diameter was 5 × 10⁻⁵ m, (Pevere et al. 2006). The properties of the lateral walls for the soil were as follows: ρ 2.650 kg·m⁻³, Cp 1,585 J·kg⁻¹·K⁻¹ (Evett et al. 2012) and λ 1.02 W·m⁻¹·K⁻¹ (Garcia 2003). The geomembrane properties were: ρ 925 kg·m⁻³, Cp 2,300 J·kg⁻¹·K⁻¹ and λ 0.33 W·m⁻¹·K⁻¹.

The governing equations were based on the Navier–Stokes equations solved in the transient state for the fluid and the dispersed phase. They were used in three dimensions (see Equations (1)–(6)). The model and its equations were addressed using a Eulerian multiphase approach because it offered acceptable results compared with the experimental data and lower computational resources than the Lagrangian approach.

In this study, the constants of the default turbulence model were assumed. These constants had the following values: C₁ɛ = 1.44, and C₂ = 1.92.

The solar load model was characterized by the parameters of the geographical location for the research station in Ginebra. The time zone was −5 GMT. The geometry was oriented according to the actual geographical coordinates as follows: Z-axis = north, X-axis = east. More details on the formulated equations can be found in the literature (Versteeg & Malalasekera 2007; ANSYS 2013).

To simulate the effect of the wind on the surface of the lagoon, the prevailing direction of the wind was used with an average speed of 0.81 m·s⁻¹ and northeast direction. Additionally, the penetration of incident solar radiation was simulated with the model and compared with experimental data measured with a Secchi disk.

The results of the tracer test in the CFD model were compared and validated with the experimental results in terms of the residence time distribution curves, following the recommendations suggested by Alvarado et al. (2013). Figure 4 shows reasonable consistency between the distribution of the curves in the experimental study and the curve obtained through the simulations of the CFD model using the chemical species model. It is consistent with the proposed interpretations of Ramesh & Nilesh (2015). Figure 4 shows that the first peak of tracer concentration obtained with the CFD model appears at the same time (2 h) as in the field study. The Shapiro–Wilk test (n > 50) showed that experimental and simulated data series did not correspond to normal distribution (P < 0.05). Then, the non-parametric Levene test was applied to compare the variances of the experimental (σ2 = 4,602) and simulated (σ2 = 4,607) series, and it was found that there were no statistically significant differences between them (P < 0.05).

With the model, the water body temperature profiles were simulated. After the fifth day of the simulation, the scenario in which the geomembrane was used showed a value close to the experimental data for the effluent temperature (corresponding to the liquid leaving the pond). The average simulated value was 25.15 °C, while the experimental value was 24.77 °C, which represents a percentage of error of 1.50%. The model in which the geomembrane was not used reached an effluent temperature of 24.45 °C, with a percentage of error of 1.29%. The results of the accuracy tests showed that errors obtained for AMD% and RMD% between the model and the experimental data for both scenarios corresponded to 2% and 11% respectively. The SSE was 0.330 and 0.270 respectively. The RMSE presented error percentages of 0.165 for the experimental data and 0.063 for the CFD model. Finally, the RMSD showed error percentages of 0.170 and 0.072 for the experimental and simulated data, indicating that the model presents less variability and greater robustness in the results. The two simulated scenarios gave results that were very close to the experimental values. However, the geomembrane allowed an effluent temperature increase of 0.7 °C with respect to the soil. This increase would be favored by the black color of the geomembrane and its greater heat capacity.

To validate the temperature profiles, three points inside the FL corresponding to point P₂ and its three depths were compared with the simulated data of the CFD model (see Figures 5–7).

Figures showed good concordance between the experimental and simulated data. The Shapiro–Wilk test (n > 50) was applied and it was found that data series did not have normality (P < 0.05). Then, the non-parametric Levene and Kruskal–Wallis tests were applied to compare the variances and the distributions of the experimental and simulated series and it was found that there were no statistically significant differences between them (P < 0.05).

However, the geomembrane favors a higher temperature at each depth. For 0.05 m and 0.45 m, the increase was 0.7 °C, and for 1.40 m, it was 0.9 °C. Although these increments appear to be non-significant, several studies have shown that small increases in the wastewater temperature favor an increase in the efficiency of bioprocesses such as methanogenesis (Spachos & Stamatis 2011). Additionally, simulations of the energy transfer through the surrounding walls of the FL show that the geomembrane mitigates heat losses by acting as thermal insulation and transporting the solar energy from the wastewater surface to the bottom. The simulations showed that there is greater heat loss when the geomembrane cover is not present. The maximum value for the heat flux with the geomembrane was 5.90 W·m⁻² while with the soil it was 16.03 W·m⁻². This demonstrates that during the design of FLs, the selection of materials for the surrounding walls is important. The proper selection allows the designer to influence the increase in the temperature of the body of water to favor the kinetics of the bioremediation processes.

The volumetric temperature profiles in Figure 7 show that the temperature increases from the inlet to the outlet (longitudinal), and decreases from the surface to the bottom (cross-sectional). In general, the variations of the temperature are small, which demonstrates the thermal stability of the body of water. The variations can be explained by the penetration of incident solar radiation into the body of water and the solid transport phenomenon. In concordance with Sukias et al. (2001) the suspended solids disperse the solar radiation and decrease its penetration within the FL, affecting the heat transport in the body of water.

In the case of the cross-sectional gradient, it is evident that at a greater depth, there will be less penetration by this radiation. This trend is due to the characteristics of wastewater (its transmittance) and because of the higher density of suspended solids. Therefore, the bottom is expected to be the coldest area of the lagoon. Figure 8 shows that the bottom (depth 1.40 m) is the coldest zone, with a difference of 0.7 °C for the geomembrane and 0.9 °C for the soil relative to the surface. The transverse variation shown by the model coincides with experimental studies that exhibit this same temperature behavior in secondary FLs (Beran & Kargi 2005; Sah et al. 2011).

Authors such as Torres et al. (1999) stated that the temperature gradient in an FL is determined by the size of the lagoon and the environmental conditions in the lagoon location area, including day–night and seasonal cycles. In that sense, Ginebra-Valle del Cauca is a municipality located in a tropical area and characterized by stability in the environmental conditions such as the temperature, the values of which are generally maintained in the range between 19 °C and 30 °C.

It indicates that these environmental variations, characteristics of the tropical regions, are not enough to generate stratification due to the changes in the wastewater density. These stratifications alter the flow patterns in the lagoon, affecting its performance and reducing the efficiency of its operation. This stratification gives rise to three zones in the body of water of the lagoon as follows: the epilimnion, thermocline, and hypolimnion.

In this research, it was evident that experimental and simulated temperature data variations between the surface and the bottom were less than 5 °C. In concordance with Torres et al. (1999) these differences can be viewed as the absence of a gradient due to the wastewater density changes. Therefore, it is possible to discard the effect of the temperature on the flow patterns by the change in fluid density and the absence of the three zones mentioned. In the case of the longitudinal gradient, it is influenced by the solid transport phenomenon. Authors such as Shilton & Harrison (2003) have shown that solids with a larger particle diameter (80 μm) tend to be deposited near the entrance of the tributary, whereas smaller ones (20–50 μm) are transported by the fluid to different areas of the lagoon such as recirculation zones or all the way to the output of the effluent. This situation generates a longitudinal concentration gradient that should allow a greater entrance of solar radiation into the final section of the lagoon. In this study, the experimental and simulated data show that there are not longitudinal gradients for the suspended solids concentrations. The absence of these gradients can be attributed to the conventional configuration of the lagoon and its length–width ratio (L × W), which is 2.56. Authors such as Abbas et al. (2006) argue that the best treatment efficiencies occur in lagoons with an L × W ratio equal to 4 with two or four baffles at L/3 (two baffles) or L/5 (four baffles).

By simulating the effect of the wind on the surface of the lagoon an increase in the temperature of the effluent and point P₂ was observed. The simulated effluent temperature was 26.65 °C, which was found to be 1.88 °C above the experimental temperature. With these data, the percentage of error in the effluent was 7.6%. It would indicate that wind on the wastewater surface destratified the liquid column and improved the mixing, affecting the hydrodynamics of the lagoon (Sweeney et al. 2007). Table 1 presents the comparison between the experimental and simulated temperature data for the geomembrane without and with wind at point P₂. It is possible to observe an increase in the temperature in the simulations with wind, with greater percentages of errors at depths of 0.45 m and 1.40 m than on the surface. At these depths, the error was 5.77% and 7.59%, respectively. These error percentages were considerably lower than the 15% error obtained by Kellner & Pires (2002); considering the several factors that occur in the thermal phenomena in waste stabilization ponds, those differences can be considered small.

Additionally, these values were significantly lower than the error of 19.05% obtained in a study performed by Ukpong (2013b). In that study, Ukpong developed a model using a code in the programming language Fortran. The differences in error percentages are an indicator of the benefits offered by CFD over other similar modeling techniques.

Wind can have an important effect on the behavior of lagoons; this is because it favors the vertical mixing of the lagoon components, improving the behavior of a fully mixed reactor (Gu & Stefan 1995). However, three factors are essential when analyzing the effect of wind on the behavior of the lagoon: the surface area, and wind speed and direction. Ukpong (2013a) showed that the effect of three distinct velocities (2.27, 1.88, and 1.64 m·s⁻¹) on a lagoon did not produce significant differences in the concentration of chemical species and temperature on the surface, but it did lead to small differences in the middle depths and at the bottom. Similar behavior can be observed in the temperature results in Table 1. The Shapiro–Wilk test (n > 50) was applied to the simulated temperature data obtained with the geomembrane with and without wind. In both cases, the data series did not correspond to normal distributions (P < 0.05). The non-parametric Mann–Whitney test was applied to compare the percentages of error. On the surface there were no statistically significant differences (P < 0.05) and at 0.45 m and 1.40 m statistical differences were found (P > 0.05). These results indicated that the behavior obtained with the CFD model was similar to results obtained by Ukpong (2013a).

The model showed that the penetration of solar radiation into the body of water reaches a 0.12 m depth, a value that is within the range from 0.11 m to 0.14 m of the penetration measured experimentally with a Secchi disk. To improve the results of the CFD model, it is recommended that different wavelengths of the solar radiation, UV-A, UV-B and PAR (photosynthetically active radiation), are evaluated and compared with experimental data as suggested by Passos et al. (2019). The penetration profiles of incident solar radiation into the body of water are shown in Figure 9. In the figure, it is possible to observe the attenuation of the solar radiation with increasing depth. This attenuation is primarily due to the absorption of light by water and humic substances Ragush et al. (2017), the dispersion due to suspended particles (Sukias et al. 2001) and the climate where the ponds are located (Heaven et al. 2005). These factors lead to longer wavelengths penetrating much better than short wavelengths (Curtis et al. 1994); for example, Dias & von Sperling (2017) reported that in a maturation pond operating in a tropical climate UV-A and UV-B radiations were detected at 10 cm from the surface, but from 15 cm both were undetectable, while PAR was still detected at 30 cm of depth. In this context, the model demonstrated the attenuation of solar radiation within the FL at 0.12 m of depth.

The comparison between the simulated data obtained with the CFD model for the concentration of suspended solids and solar radiation penetration in the wastewater at points P₁, P₂ and P₃ and their respective depths demonstrated that higher penetration values of solar radiation coincide with solid concentrations less than 0.026% w/v. Therefore, design improvements of FLs should serve to improve sedimentation of suspended solids and facilitate the penetration of solar radiation into wastewater.

A new method of including parameters such as temperature profiles, heat transfer, solar radiation penetration and the effect of wind conditions, and ways of measuring them, to characterize the lagoons and simulate their behavior in a computerized environment was developed.

The validated CFD model could be used as a support tool to asses design improvements with the intent to increase the energy performance of FLs. The model was used to assess two types of material for the surrounding walls showing that the proper selection allows the designer to influence the increase in the temperature of the body of water to favor the kinetics of the bioremediation processes.

The authors thank the staff of the Research Station of Wastewater and Reuse in the municipality of Ginebra-Valle del Cauca, Colombia, for the support they provided.

The Universidad del Valle and Colciencias and their credit-grant program, administered by Colfuturo, supported this work.