Detailed experimental and numerical analysis of hydrodynamics in the outflow measurement channel of a sewage treatment plant

Water is an invaluable resource, essential for all living beings. The continuous increase of the population in the cities has led to the degradation of many natural ecosystems such as rivers and oceans, partly as a consequence of the discharge of improperly treated urban and industrial wastewater. Therefore, from a health and environmental point of view, it is of utmost importance to correctly ensure the sanitation of urban centers and the purification of their wastewater through sanitary works.

The requirement that the treated discharge meets quality standards in accordance with the assimilative capacity of the receiver water body drives the development of new technologies and optimization of the existing treatment systems in order to improve performance in the elimination of pollutants and the sustainability of treatment systems.

Commonly, hydrodynamic singularities are often observed in the turbulent flow present in water and sewage treatment systems that can affect the plant's operation and efficiency and in consequence the final quality of the treated effluent. This applies to all treatment processes of a plant (primary treatment, settling, chlorination, etc.). Therefore, a complete characterization of these singularities is required and for this, it is necessary to use experimental and numerical techniques with appropriate spatial and temporal resolutions. Regarding the experimental techniques, during the last decade, the use of an Acoustic Doppler technology started to be relevant for experimental in situ measurements in water and sewage treatment systems monitoring and characterizing the flows. Garcia & Garcia (2006) have characterized the flow turbulence and the induced mixing by air bubbles–plumes in sewage tanks under non-stratified conditions using an Acoustic Doppler Velocimetry (ADV). Kiss & Patziger (2013) used an ADV to better understand the flow patterns in the primary sedimentation tank of a wastewater treatment plant (WWTP). In addition, an Acoustic Doppler Current Profiler (ADCP) has been used to characterize flow in clarifiers and key physical processes: hydrodynamics, particle aggregation, and floc breakup (De Clercq et al. 2002; Vanrolleghem et al. 2006). In addition, during recent years, the implementation of Large-Scale Particle Tracking Velocimetry (LS-PTV) has experienced an important development. This technique has been successfully used for the characterization of turbulent flows in large experimental laboratory facilities (for example, in physical models: Patalano & García 2016, on river bends: Patalano et al. 2014) and in a first approach to characterize the incoming turbulent flow to the treatment plant decanters and to validate a Reynolds-Averaged Navier–Stokes Equations (RANS) model to study alternatives for operation (Ragessi 2019).

On the other hand, computational fluid dynamics (CFD) has become an alternative for the study of turbulent flow in treatments units of water and wastewater treatment processes (Goula et al. 2008; Shahrokhi et al. 2011; Zhang 2014; Patziger & Kiss 2015; Wei et al. 2016; Ragessi et al. 2019). In the last 5 years, the amount of available references regarding good practices or protocols in the use of CFD models has grown significantly, applied not only to improve the design but also to daily monitoring and operation of treatment units. In spite of this, the state of the art indicates that the characterization of hydrodynamics, turbulence, flow instabilities, and their effects on components of treatment plants is still an engineering challenge.

In this way, Laurent et al. (2014) present a protocol for using CFD as a practical support tool for everyday use in WWTPs. Also, Patziger et al. (2016) reported the results of an investigation that combines experimental techniques and CFD, as part of a whole research program for the design and operation of primary settlers of WWTPs.

Samstag et al. (2016) provide a general vision of CFD application in the treatment units of water and wastewater. In their document, they articulate the state of the art and vacancy areas for future research. Furthermore, they conclude that in general (except possibly in the disinfection process) numerical models have not been used in a generalized or routine way as a design, risk management, or problem-solving tool. Finally, Wicklein et al. (2016), as part of the research work of the International Water Association (IWA), have synthesized the main characteristics of good and modern numerical modeling practices applied to wastewater treatment.

The aforementioned study indicates that the flow turbulence characterization in water and WWTP components is still a challenge and that it is of great convenience to complement modern experimental techniques with numerical tools in order to obtain a more detailed flow characterization.

For this reason, this paper presents a study in the direction of this challenge. In particular, three-dimensional flow patterns (i.e. flow instabilities and recirculation zones) with high spatial and temporal variability are studied in the outflow measurement channel of a sewage treatment plant which complicates the accurate quantification of the outflow discharge. The precise data of the flows present in the plant are very important to optimize its operation. Therefore, it is necessary to optimize the design of these measurement channels using advanced experimental and numerical techniques suitable for this type of study. Combining optimal experimental and numerical techniques is a challenge.

This paper presents the case study of the ‘Bajo Grande’ Wastewater Treatment Plant (BG-WWTP) which is located at the eastern end of the city of Córdoba, Argentina - one of the cities with the highest population growth in the country in recent years. The treatment plan discharges the effluent into the Suquia River whose annual average flow rate is 10 m³/s with a clear seasonality of the drained flows with significantly higher values in the rainy period (October–April) but with values of the order of 1 m³/s in the dry period. The design capacity of the BG-WWTP is 10,000 m³/h (2.78 m³/s), which clearly shows the importance of the effluent of BG-WWTP because any deficiency or problem in the treatment could have a high impact on the water quality of the river. It should be noted that the Suquia River flows into a lake system (named ‘Laguna del Plata’) which is connected to the Mar Chiquita lake in Cordoba province. This lake is currently protected by national and international regulations, due to its important biodiversity.

This paper explores the flow conditions in the outflow measurement channel of the BG-WWTP downstream of the contact chamber (prior to its discharge into the Suquia River) as indicated in Figure 1(a). The outflow measurement system, which assumes a steady and homogenous behavior of the flow, is installed in this channel. However, at the exit of the contact chamber, a strong recirculation pattern was observed affecting the whole outflow measurement channel, due to the presence of an abrupt flow separation induced by a 180° sharp bend at the inlet section of the channel (Figure 1(b)). The shear layer that limits the recirculation zone presents time instabilities that are manifested as low-frequency fluctuations in the local flow velocities in the flume where the flow measurement system is installed, which is not recommended because the effective flow area in the cross-section (area with positive flow velocities) changes over time and the measured flow rates with conventional methods have great associated uncertainty. This uncertainty in the outflow discharge makes the optimization of the biological and chemical treatment difficult and, besides, the quantification of the potential effects on the receiving water body, which is the main objective of the treatment system.

In this work, first, two experimental techniques were used for detailed characterization of the turbulent flow in the flume: (1) an ADCP that records vertical profiles of the three velocity vector components in the water column when it is stationary and the flow discharge when it is installed in a moving platform and moved in the cross-section and (2) the LS-PTV, a non-intrusive approach to measure instantaneous velocities at the free surface of a water body. The application of these techniques (which were not initially designed to be used in treatment plants) allowed us to record flow information with high spatial–temporal resolution. The detailed hydrodynamic characterization allows a better understanding of the effects of turbulence and flow instabilities in the study zone and also provides a basis to validate the numerical model.

In the second stage, a three-dimensional numerical model was used. The Detached Eddy Simulation (DES) technique, which combines the Large Eddy Simulation (LES) technique with the RANS technique, was used to model the turbulent flow. This strategy is convenient in flows with high Reynolds number, flow separation, and recirculation zones that change over time. According to Koken et al. (2013), for most complex turbulent flows in which strong adverse pressure gradients are present (as observed in the U-shaped curve), DES has been shown to be more accurate compared not only to RANS and unsteady RANS but also to LES with wall functions (Spalart 2009; Constantinescu et al. 2011).

In this case, the free open-source code OpenFOAM^® (Open CFD) was used and the details of the numerical model are described in the DES model section.

In the present study, the DES approach was preferred (and not the well-resolved LES) because the main goal of the present paper is to perform Eddy-resolving simulations at Reynolds numbers corresponding to field conditions. The use of hybrid RANS-LES models, like DES, is presently the best option to accurately simulate flow both at laboratory- and field-scale Reynolds numbers, with an acceptable computational cost.

This approach makes it possible (1) to assess the predictive capabilities of DES to capture the details of the mean flow developing over the entire outflow channel and (2) to use the numerical results to examine details of the three-dimensional flow (e.g. the time evolution of the recirculation zone width) and turbulence structure and estimate relevant quantities (e.g. turbulent scale, power spectrum in frequency domain) that can be difficult to estimate with in situ experimental techniques.

A total of three field campaigns for three different operating conditions of the plant were performed in the outflow channel indicated in Figure 1(b) of the BG-WWTP. The outflow channel has a rectangular cross-section with a width of 6.5 m and is lined with concrete. The main goals of these measurements in this zone were as follows: (a) quantifying BG-WWTP outflow discharges and (b) characterizing the temporal and spatial evolution of the flow velocity fluctuations (describing the flow instabilities and the recirculation zones).

The measurements were made using an ADCP and the LS-PTV technique. Both techniques have the advantage of being able to characterize the flow field with high spatial–temporal resolution. An ADCP records one vertical velocity profile per second (recording frequency 1 Hz) and finally integrates the velocity field throughout the cross-section to calculate the flow rate. Instead, LS-PTV spatial and temporal resolution depends on the digital camera, but generally speaking, it is HD (1,280 px × 720 px) and 30 fps or higher. Thus, once validated, LS-PTV has the advantage of being a non-intrusive technique to measure the instant cross-section velocity profile and the mean cross-section velocity profile. This is a relative advantage over other devices that measure the water level and velocity in a single point (commonly used in channels) or mechanical current meter, since LS-PTV calculates the flow rate from the surface velocity profile of the entire cross-section.

A disadvantage of using LS-PTV is that the real discharge must be calculated as a function of the discharge estimated with the surface velocity, for which a proportionality coefficient must be adjusted. This may not be an important issue for simple flows, but this is not the case. For this reason, the ADCP data were used to validate the LS-PTV results on the outflow channel (see LS-PTV measurements and Outflow discharge measurements sections).

The three operating conditions correspond to a medium flow (below 50% of the treatment plant capacity), an intermediate flow, and a flow close to the design capacity of the plant. Smaller flow discharges were not characterized because it is not possible to survey these conditions due to the continuous demand of the sewerage network.

Figure 1(b) indicates the discharge-measurement section located immediately downstream of a pedestrian walkway, 34 m downstream from the 180° sharp bend vertex. An ADCP ‘RiverSurveyour S5’ mounted on a mobile platform was used for these flow discharge measurements. Based on flow depth and velocity, and flow turbulence levels, the software ‘RiverSurveyour Live’ (RSL) optimizes the recoding configuration and adapts the acoustic pulse in order to obtain the highest spatial resolution possible: with cells up to 2 cm.

In all the measurements, the Bottom Track channel in RSL was used as a reference system (Bottom Track in RSL). A potential velocity law was adjusted over the entire velocity profile (exponent=1/6) to estimate the flows in the unmeasured regions: (1) near the free surface of the flow (due to the needed 7 cm submergence of the instrument) and (2) near the bottom of the channel (to avoid considering effects of the bottom in the acoustic pulse). On the other hand, the unmeasured flows in the right and left banks were estimated using the recommended method by RSL for vertical riverbanks.

Flow instabilities along the flume made it difficult to determine the flow discharge (from a mobile platform). The number of transects and the measurement exposure time (duration of the measurement) are critical factors in reducing discharge-measurement uncertainty and depend on the flow conditions (flow depth, flow velocity). Therefore, in order to minimize the uncertainty in the measured values, multiple transects were performed instead of the four generally recommended to estimate flow discharge (Tarrab et al. 2012; Rahimpour et al. 2021). The number of performed transects exceeds the minimum exposure time required for ADCP measurements (12 min, Mueller et al. 2013).

Finally, the flow discharge was calculated as the average flow of the multiple transects: 20 transects, 12 transects, and 16 transects for in situ measurements C1, C2, and C3, respectively, as indicated in Table 1.

Table 1 shows the main characteristics of the experimental conditions in Section A-A of the flume (Figure 1(b)) and the sampling characteristics of the ADCP measurements. In the table, H indicates the depth of flow, indicates the recording frequency of the ADCP, indicates the number of cells in the vertical, and P indicates the number of profiles in the transect.

The ADCP register frequency is fixed for the selected model, and the cell size is automatically adjusted with de Riversurveyor software to optimize performance and resolution. The average maximum number of ensembles per transect (⁠⁠) depends on a tethered boat velocity and the average boat speed during each transect normally should be less than or equal to the average water speed.

The analysis of the ADCP data was performed with the Qrev program (Mueller 2016).

The sampled flow discharges were compared with values indirectly estimated as a function of the operating frequency of the pumps at the pumping station (as it was used for the company in charge of the operation of the BG-WWTP based on the issues observed in the outflow measurement channel).

In addition to the discharge measurements performed using ADCP, the LS-PTV technique was used in the channel area immediately upstream of Section A-A (highlighted rectangle in Figure 1(b)) to characterize the flow field. A camera with a video resolution of 1,280 px × 720 px and a 30 Hz image capture rate (30 fps) was used. A video of 90-s duration was analyzed.

In the implementation of the technique, it was necessary to rectify the images (Figure 2(a)) because the camera was not orthonormal to the flow field. For the mentioned rectification, four control points were marked on the channel banks. See the rectified image of the recording area in Figure 2(b).

Laboratory tests were performed to determine the optimum tracers used for LS-PTV technique implementations (Patalano et al. 2017). After the tests, elements formed by two plastic circles adhered to each other of 65 mm diameter were selected (see Figure 2(c)). The reasons for this choice were mainly its visibility and low cost. The software used for particle tracking, PTVlab, is developed by coauthors of this manuscript and has the particularity to track white and circular particles. The effects on the measurements of the type of tracers used have been evaluated by comparing them with other tracers.

The tracers were seeded manually from upstream to the measurements section to achieve a homogeneous distribution on the surface of the flow. In addition, the configuration achieved a weight such that the action of the wind does not modify its speed and at the same time maintains buoyancy and constitutes symmetrical particles with respect to its two faces. The tracers weighed 3 g, and on the day of the measurement, the wind speed was less than 15 km/h at 10 m height (recorded by the international airport 20 km away from the plant). On the day of the experiment, the wind was barely perceptible without any gusts and the authors decided to neglect its drag effect on the tracer.

The flow rate was estimated from the mean surface velocity profile and the surveyed section (width and depth).

The ratio depends on many parameters, such as the geometry section, the bed roughness, the geometry upstream and downstream of the section, secondary flows, and others. It is up to the user to choose the most reliable value of ⁠. Cheng et al. (2004) concluded that it is reliable to use the surface velocity as an index to determine a river flow discharge because the value of always falls in the same range from 0.80 to 0.93.

In this case, α was calculated based on the vertical velocity profile fitted from all ADCP data processed with the Qrev program. In addition, α was also calculated using the velocity field solved with the DES model. Then, the flow rates calculated with ADCP and LS-PTV (using both α values) are compared in the Results and Discussion section (see Outflow discharge measurements section).

Comparisons were performed of the velocity flow field sampled in multiple transects: 20 transects, 12 transects, and 16 transects for in situ measurements C1, C2, and C3, respectively, as indicated in Table 1.

In addition, to characterize the temporal evolution of the flow velocity and instabilities, during in situ measurements C2 and C3, velocity profiles were recorded for 10 min with the same ADCP on a fixed platform located at Point B: center of the Section A-A (Figure 1(b)). The times series of vertical mean longitudinal velocity profiles were processed to determine the time scale of the fluctuations by analyzing the power spectrum in the frequency domain, and autocorrelation function, among others following a similar methodology presented by Garcia & Garcia (2006).

The LS-PTV technique was used to characterize the cross-sectional velocity profiles in the channel area immediately upstream of Section A-A (the same section where ADCP measurements were performed).

The experimental information was only able to be obtained in one area of the channel because it was only possible to use the ADCP and apply LS-PTV in a zone of the outflow channel. Hence, in addition to in situ measurements, a CFD simulation of the entire channel was made for the purpose of capturing the flow separation and the instabilities that are manifested as low-frequency fluctuations in the local flow velocities.

A three-dimensional numerical model was used, which combines the LES and RANS Equations, known as DES: the unsteady RANS models are applied in the boundary layer, while the LES model is used for the separated regions.

The DES model used in this work is based on the one-equation modified turbulence viscosity (⁠⁠) SpalartAllmarasDES model for incompressible flow, with the low Reynolds number term correction function in the proximity of the walls.

is a laminar suppression term and fully turbulent flows at high Reynolds numbers takes small values.

The default value of the DES constant was calibrated by OpenFOAM^® (OpenCFD) under the hypothesis of isotropic turbulence and the results agree with the suggested values in the referenced publications (Spalart et al. 1997; Spalart et al. 2006). The default values are present in Table 2.

The DES approach, like the one used here, has been used to represent and study the hydrodynamics of flow in curves (such as the flow observed in the BG-WWTP flume). In particular, by referring to curved channels Van Balen et al. (2010), Constantinescu et al. (2011) and Koken et al. (2013) have managed to model three-dimensional flow separation in the inner margin with satisfactory results.

The temporal PIMPLE scheme (pimpleFoam in OpenFOAM^®) was used, which merges the PISO (Issa 1986) schemes and SIMPLE (Patankar 1980) for non-stationary flows, for constant density and viscosity. In general, better stability is obtained from PIMPLE over the PISO algorithm with the use of outer correctors, especially when dealing with large time steps where the maximum Courant number may consistently be above 1 or when the nature of the solution is inherently unstable (and needs a long-time simulation). In this case, one outer-corrector and three inner correctors (to correct pressure in each iteration) were set up mimicking a PISO scheme, with a limit for the max Courant number of 1. The time step was approximately 1×10⁻³ s and the results were saved (writing control) every 0.04 s (25 Hz frequency).

The use of freely-available open-source codes, like OpenFOAM^®, allows a continuous community-based improvement of the model and avoids having to pay for costly software licenses. In this sense, openFOAM(R) has been widely used and validated in multiple areas of engineering and research, both by industries and universities (Jasak et al. 2007; Jasak 2009; Lysenko et al. 2013; Verstraete et al. 2013).

In the detailed simulation, a mesh of variable-size cells was used. Away from the zone of separation and recirculation of the flow, volumes of 4 cm of the side were used, and in the area with higher velocity gradients, near the walls, cells of 0.5 mm of the side were used. The mesh is composed of 11.226.460 cells. In Figure 3, the domain and the detail of the resolution in the area near the 180° sharp bend are shown.

Several approaches are presented in the literature to evaluate the resolution of LES (or DES). In general, it is recommended that in a good LES simulation, the integral length scale L should cover 8–16 cells. Other authors, like Pope (2004), suggest that when 80% of the turbulent kinetic energy is resolved, the LES can be considered to be well-resolved. Davidson (2009) and, more recently, Seifollahi Moghadam et al. (2021) reported that the ratio of the resolved turbulent kinetic energy to the total one and the energy spectra are not a good measure of LES resolution in different test cases.

Kuczaj & Komen (2010) have shown that mesh-size estimation based on the Taylor micro-scale a posteriori computed from LES (⁠⁠) gives sufficiently accurate results. Addad et al. (2008) and Kuczaj et al. (2010) used a precursor RANS model in order to examine characteristic turbulence scales that can be estimated and used for LES mesh construction. They conclude that meeting this requirement allows estimating turbulent flow characteristics with sufficient accuracy.

They used for the minimal LES mesh resolution (recommended by Addad et al. 2008).

In this case, a precursor RANS model resulted in that: and ⁠. The largest side of a cell in the used mesh is between this range (4 cm), while near the curve and the walls the resolution is 0.5 mm (⁠⁠).

More recently, Seifollahi Moghadam et al. (2021) indicate that the two-point correlations are suggested as the best measures for estimating LES resolution. But at least the results of three meshes should be compared in a mesh sensibility analysis (e.g. Root Mean Square (RMS) of velocity fluctuations). This approach would have a relatively high computational cost and would be more time-consuming for a prototype-scale model. In future works, this approximation will be used to evaluate the impact on the estimation of turbulent flow properties due to the increase in the resolution of the mesh.

Gmsh (Geuzaine & Remacle 2009) was used to resolve the mesh. Gmsh is an open-source 3D finite element mesh generator with a built-in CAD: the software provides a fast, light, and user-friendly meshing tool with parametric input and advanced visualization capabilities. These characteristics, added to the fact that the software is open-source and compatible with openFOAM(R), are the reasons for selecting Gmsh to solve the mesh. To import the mesh the gmshToFoam application was used.

Table 3 summarizes the boundary conditions imposed in the computational domain for the velocity vector (U), the pressure (p), the turbulent kinetic energy (k), and the turbulent Kinematic viscosity (⁠⁠, nuTilda en OpenFOAM^®).

In the prototype, the free surface did not show variations in level, and in the numerical model, a rigid surface was imposed without shear stress so that it did not affect the internal flow (slip conditions). In this way, it is possible to simplify the two-phase modeling scheme (water–air) by simulating the free surface of another of the (water) phases. This simplification allows us to reduce the necessary computational effort and work with a less complex modeling scheme with respect to the boundary conditions, physical parameters, and equations to solve.

The validity of the rigid-lid hypothesis was then verified using the pressure field values at the top surface of the model. It was observed that the maximum relative difference does not exceed 10% in the bends area.

At the inlet, a constant uniform value for U was imposed to obtain a constant flow rate of 2.22 m³/s. It was difficult to generate turbulent inflow fluctuations that are realistic for this specific problem, and using a constant velocity value at the inlet patch is a simplification of the real flow inlet condition. Thus, to lessen the effect of using on the model results, the entry edge condition was moved away and a first serpentine tank U-bend was incorporated. This allowed the development of turbulence downstream of the first U-bend as a consequence of the shear layer due to the flow separation. Other works have shown that the inlet boundary condition has apparently no influence on the flow field downstream of an obstacle and in the separation zone and that in this area the mesh resolution is more relevant (Montorfano et al. 2013).

The outlet boundary was controlled by a zero gradient (Neuman BC) velocity condition and no-slip boundary conditions were imposed on the side and bottom walls. The outlet path was placed 3 m upstream of the outlet ducts at the end of the outflow channel: when the flow changes its regime from subcritical to supercritical (see Figure 1(c)).

The standard wall function for k was employed so that the governing equations could be solved correctly. The turbulence model wall functions have restrictions on the non-dimensional distance (⁠⁠) value at the wall. In this case, it was verified that the cells have the proper size near domain walls, checking that values were between 30 and 100 and thus it is valid to use this wall function.

Table 4 shows the flow rates measured with ADCP (⁠⁠), the discharges computed as a function of the frequency of the pumps (⁠⁠), and the relative error (⁠⁠) based on the values measured with ADCP.

In general, the use of the information reported by the pumping system during all three in situ measurements overestimated the discharged effluent. However, Table 4 shows good agreement between the estimated values during in situ measurement C3 by both monitoring systems (pumping system and ADCP). The relative differences observed in this in situ measurement were less than 3%. A larger difference is observed for the lower flow rate: the relative difference is 67 and 19% for the in situ measurements C2 and C3, respectively. The outflow discharge estimated by the pumping system presents an inconsistency when the corresponding flow depths in the flume are analyzed: of the in situ measurement, C2 is greater than C3 but has a lower flow depth (⁠⁠). According to ADCP data, the treatment plant was operating at 44, 70, and 80% of its full capacity (2.78 m³/s) during in situ measurements C1, C2, and C3, respectively.

To calculate the flow rate with the LS-PTV technique the value was estimated. First, was calculated based on the vertical velocity profile fitted from all ADCP data processed with the Qrev program, and, second, the value was adjusted based on CFD model results, using the ratio between the mean flow velocity and the surface velocity. The ratio was equal to 0.95 and 0.97 from ADCP data and DES model results, respectively.

The flow rate calculated by the LS-PTV velocity profile turned out to be 2.27 and 2.32 m³/s (using both values). The relative difference between ADCP flow rates in C3 (2.22 m³/s) is 2.5 and 4%, in both cases, less than 10%.

The ratio is slightly higher but close to recommended values (0.8–0.93). But the results show that the LS-PTV technique can be considered a cost-effective alternative and simple application to carry out local flow rate measurements.

Using the ADCP flow velocity measurements from the moving platforms, recirculation zone width (⁠ see Table 4) was measured from the left margin of section A to the point where the zero velocity occurs (changes from negative to positive flow). The recirculation zone affects from 20 to 30% of the cross-section width: the recirculation area extends to 2, 1.8, and 1.6 m from the left bank in C1, C2, and C3, respectively (see Table 4). This would indicate that for lower flow rates than the WWTP-BG design capacity, the instabilities effect is larger. While for flow rates of the order of the plant design condition, velocities increase, inertial forces gain relative relevance as is shown by the Froude and Reynolds numbers.

The Froude number (⁠⁠) in the flume was very low in the three in situ measurements (⁠⁠). In this case, u is the longitudinal velocity component, and it is calculated as a function of the measured flow rate and the area of the flume cross-section (⁠⁠), g gravity acceleration, and H height of water surface measured (indicated in Table 4).

The Reynolds number was calculated as ⁠, where the hydraulic radius (⁠⁠) is the ratio between the cross-sectional area of the channel divided by the wetted perimeter and is the kinematic viscosity equal to 1×10⁻⁶ m²/s.

The low numbers facilitate the propagation of instabilities along the flume and do not allow the development of a uniform flow. In this type of complex flow, it is difficult to quantify the discharge. The development of recirculation zones and shortcuts along the chamber would indicate that there is a wide range of contact times distribution, that must be taken into account to adjust the chlorine dosage in the disinfection treatment.

An example of this can be seen in Figure 4(a) where the longitudinal velocity component (⁠⁠) was plotted for Section A-A. The color plot is the result of the averaged transects for in situ measurement C3. On the right bank, the flow velocity is of the order of 1 m/s, while on the left margin it reaches negative values (which indicates recirculation) of the order of 0.5 m/s in the entire flow depth. The highlighted area inside the dashed line indicates the area where the velocity changes from negative to positive values (shear layer zone). Results for in situ measurements C1 and C2 report a similar flow behavior.

In Figure 4(b), DES model results were plotted. The maximum and minimum values of are similar to ADCP data (1 and–0.5 m/s), and the recirculation zone affects 30% of the cross-section width. These results are very similar to those reported by the averaged transects registered with ADCP in C3.

Figure 4(c) shows the relative difference between ADCP and DES results for the longitudinal velocity component values. The major differences are a little higher than 10% near the shear layer, and from 2 to 6% in the recirculation zone. These differences are not due to changes in the flow discharge over time.

The shear layer is a very dynamic zone, where recirculation width changes over time (therefore the maximum and minimum flow velocity values also change) due to large time scales with periods ranging from 100 to 10 s in the macro-turbulence process (this will be present below in detail, see Figure 7). This can explain the relative difference between ADCP and DES results for the longitudinal velocity component values, but so far it can be said that the DES model approximates mean flow very well.

Another way to view this change in recirculation zone width is through the instantaneous velocity values across the section. Figure 5(a) compares ADCP velocity data with DES model results at a distance from the water surface of 60% of flow depth in Section A-A. In this figure, the blue points are the instantaneous velocities values from transects performed with ADCP in C3, and the red line is the mean-velocity profile (ADCP-Mean). While the gray area delimited by the dashed line represents the 95% confidence interval of the mean-velocity profile calculated with the percentile method from 10 min of DES simulation (DES-Data CI).

It can be seen that most of the instantaneous velocity values (blue points – ADCP) are within the gray area (DES-Data CI). On the right bank, low velocity values of 0.1–0.2 m/s were observed with the ADCP and the DES model. But the DES model reported higher velocities than ADCP: 1.4 and 1.2 m/s, respectively.

On the left bank, the recirculation zone, ADCP, and the DES model reported velocities between −0.7 and −0.6 m/s. Also, both ADCP and DES models report similar positive velocities: 0.3–0.2 m/s.

The differences may be since the results from the DES model are continuous in time and ADCP only collects instantaneous velocity data during the transects. Also, in each transect, the ADCP does not carry out exactly the same trajectory near sections A-A.

Similar results can be observed in Figure 5(b), where surface velocity data in Section A-A are plotted. In this figure, the average velocity profile determined with the LS-PTV technique was added (black line in Figure 5(b)). The recirculation zone extends up to progressive 1.0 m and maximum velocities are 0.9 m/s. Differences between mean-velocity profiles may be due to different times over which each was averaged: the LS-PTV profile is from a 90 s measurement and the DES profile is from 600 s of simulation. LS-PTV profiles are within the confidence intervals of velocity profile from DES model results (DES-Data CI), therefore, they can be considered statistically similar.

In this case, the channel has a uniform rectangular section but the measured flow is downstream of a curve. Therefore, a value for the ratio between the mean flow velocity and the surface velocity was adjusted based on ADCP results as explained above (in this case, was equal to 0.95).

In Figure 5(a), the instantaneous velocity values (blue points – ADCP) show that the recirculation zone width (and the shear layer location) changes over time: it extends up to the progressive 3.2 m (maximum reported) and can be retracted back to progressives less than 1.0 m.

A very similar flow behavior can be seen in Figure 6(a), where u longitudinal velocities at a depth from the water surface of 60% of the depth of the flow in Section A-A were plotted over time. In this plot, a clear and sharp transverse variation in velocity can be observed. The blue line indicates the shear layer position. The shear layer lies between progressives 0.9 and 2.90 m and moves between these values in time. The location of this shear layer coincides with the location of the recirculation zone reported by ADCP data in Figure 5(a).

Time-fluctuations affect the entire vertical velocity profile. This behavior was observed with ADCP data in C1, C2, and C3. This also can be seen in the DES model, not only for mean flow (Figure 4(b)) but also for the instantaneous velocity fields in Figure 6(a) and 6(b).

In addition, time–velocity signals from the DES model and a stationary ADCP measurement were analyzed.

The ADCP was placed in the middle of sections A-A (near the Point B location) to record vertical velocity profiles for 10 min with a sampling frequency of 1 Hz. Each velocity component (x, y, and z) was calculated from beams velocity data. From the DES model, velocity signals were sampled with a frequency of 25 Hz in three points near the Point B location (25 cm close to each other). The time-lapse (10 min) is enough to characterize the low-frequency fluctuations and flow temporal evolution.

Figure 7 shows the computed power spectrum in the frequency domain (G) from magnitude velocity signals, for both stationary ADCP and the DES model. The power spectrum was computed through the fast Fourier transform (FFT) of the windowed data.

The macro-scale turbulence process is characterized as shown in Figure 7 through peaks at the low-frequency component of the power spectrum, both observed (ADCP) and simulated: the larger structures with frequencies between 0.009 and 0.025 Hz (about 110–40 s) and a second peak in the bandwidth of 0.06–0.07 Hz (about 16 s). These periods coincide with the temporal evolution of the shear layer location illustrated in Figure 6.

To characterize the larger flow turbulent structures, a similar methodology used by Garcia & Garcia (2006) was used.

For the transverse component, the autocorrelation function is and for the vertical component ⁠. Here ⁠, ⁠, and are the mean values of the time–velocity signals ⁠, ⁠, and ⁠.

The concept of convective velocity is widely used when it is desired to estimate length scales of processes based on time series using the Taylor frozen turbulence hypothesis. There are different ways to estimate that convective velocity (⁠⁠). The simplest way is to use the mean longitudinal velocity (⁠⁠). However, this is not suitable for three-dimensional processes. For this reason, the methodology used in this paper follows Garcia & Garcia (2006).

Table 5 shows computed turbulence large scales that are present in the flume. Wu & Patterson (1989) suggested the use of a resultant integral length scale that takes into account the three dimensionalities and anisotropy of the flow given as ⁠. In Table 5, ⁠, ⁠, ⁠, and are presented.

The integral length scale from both the DES model and ADCP shows similar results in the order of 70–75% of the section width (⁠⁠) and 74% of the water depth (⁠⁠).

This is also verified in the behavior of transverse velocities in different sections of the DES model. Figure 8 shows the magnitude of the dimensionless transverse velocity with a mean value of the longitudinal component (⁠⁠) and the velocity vectors in four sections along the flume: placed at 5, 10, 20, and 35 m downstream from the 180° sharp bend vertex. The latter corresponds to Section A-A (Figure 1(b)).

At the exit of the curve (Figure 8(a)) the formation of a horizontal–axis secondary flow cell, with strong transverse velocity components, can be observed. The turn is clockwise and occupies 40% of the section width, on the right bank. In the next section, 5 m downstream (Figure 8(b)), the secondary flow cell is still present but the transverse components are much smaller in comparison to the mean longitudinal velocity. At 20 and 35 m downstream, in Figure 8(c) and 8(d), the clear development of the helicoidal flow is not observed but, in some cases, the transverse velocity components are of the order of 20% of the mean longitudinal velocity component. These are average results and instantaneous fluctuations may be larger, which suggests that despite the length of the flume (about 50 m) it is not possible to develop a uniform flow, and strong transversal components of velocity are observed in all its length that fluctuate in time with a wide range of scales according to the results of the spectral analysis.

This work was limited to describing in detail the turbulent flow in the outflow channel of a contact chamber (serpentine tank). The results shown here provide a basis to validate numerical models that could be used to improve the general design or to find a simple geometry for the upgrading of existing tanks with the aim of reducing recirculation zones and improved disinfectant mixing with the flow. This approach follows the same path as recent efforts to study treatment unit performance with simulation using the Reynolds-averaged governing equations in prototype-scale (Kizilaslan et al. 2019, 2020; Zhang et al. 2019) and LES model for laboratory scale (Bruno et al. 2021).

In this work, a detailed hydrodynamic characterization of the turbulent flow in the outlet flow in a WWTP was performed by means of in situ experimental techniques and a detailed numerical model. Based on the results achieved, recommendations can be made to overcome the problem of the discharge-measurement system that is operating.

First, the ADCP and the LS-PTV techniques provided essential information: quantification of the plant operating flow, flow velocity fields, the temporal evolution of the cross-section velocity profiles, and instantaneous and mean velocities fields at the free surface. Additionally, the LS-PTV technique, of low implementation cost and calibration, proved to be an alternative to performing flow rate measurements with low uncertainty.

The alternative of the LS-PTV technique for flow rate measurement has the advantage that it would not be necessary to modify the current flume geometry of the BG-WWTP and it is easy to apply.

Second, it is important to remark that the DES model was able to reproduce not only the mean behavior of the flow but also the temporal evolution of the instabilities in the outflow channel: velocity fluctuations, recirculation zone, and shear layer locations. On the other hand, the DES model proved to be a valuable tool to characterize in detail the flow hydrodynamics along the flume. Therefore, if the alternative would be to improve the inlet section geometry of the flume, the numerical model could be used in the re-design process with the objective of reducing low-frequency fluctuations and achieving the development of a uniform velocity profile: a necessary condition for discharge measurements.

Also, the larger flow turbulent structures present in the flume were characterized by time–velocity signals. From both the DES model and ADCP results, the integral length scale (⁠⁠) is in the order of 70–75% of the section width. The larger flow turbulent structures present in the flume are superimposed on the low-frequency macro-scale turbulence process with a timescale of approximately 110–40 and 16 s. These periods coincide with the temporal evolution of the shear layer location.

Third, the use of a numerical model allowed us to extend the characterization of the flow throughout the full length of the flume. The results indicate that despite the great length of the flume (about 50 m) it is not possible to develop a uniform flow, and strong transverse velocity components that fluctuate over time are observed throughout its length.

This work was limited to describing in detail the turbulent flow in the outflow channel of a contact chamber (serpentine tank). But the above indicate that DES is a valuable tool not only to characterize in detail the hydrodynamic behavior of the flow but also to incorporate it into the design process of new plants and for the upgrading of existing ones.

The numerical model could be used to improve the hydraulic design of the contact chamber that considers minimizing recirculation zones and the flow shortcuts in serpentine flow to guarantee a sufficient and uniform contact time to deactivate the pathogens. This approach will be used in future works.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.