Proposal of a new approach to perform advanced flow modeling for subsurface flow treatment wetlands

Treatment wetland (TW) is a wastewater treatment technique based on engineered systems designed to optimize processes that can occur naturally in environments and are therefore considered environmentally friendly and sustainable options for wastewater (Dotro et al. 2017). Models for wastewater treatment are often more focused on biokinetic and process-based models than in hydrodynamics of TW (Reed et al. 1995; Vymazal 2018). However, efficiency in pollutant removal in TWs is related to the system's hydrodynamics. Indeed, inadequate flow patterns can lead to operational problems in TWs, such as dead zones and preferential flow paths.

Computational fluid dynamics (CFD) has been an important tool for studying internal dynamics in TWs. The most recent models adopted the Darcy or Darcy–Forchheimer equations for modeling the porous media, which presents important details for studying TWs, such as the hydrodynamics involved inside the system (Samsó & Garcia 2013); Yuan et al. 2020; Meyer et al. 2015).

The discrete element method (DEM), on the other hand, enables the simulation of porous media at the particle scale by incorporating a vast array of particles of varying sizes and shapes, thus making it a more accurate approach for modeling granular media (Tong et al. 2020). When coupled with CFD, DEM has the potential to better represent the hydraulic behavior of TWs, making it capable of simulating flow anomalies and aspects that other models are not able to simulate.

So, in computational fluid dynamics and discrete element method (CFD-DEM) coupling, the CFD component simulates the flow of a continuous fluid based on local averaging. The DEM component, in turn, describes the behavior of granular media at the particle scale using Newton's laws taking into account the forces acting on the particles by the fluid. So, CFD-DEM is considered a robust tool especially useful for studying systems involving fluid flows and discrete particles (Golshan et al. 2020), which makes this method a promising approach for advanced TW simulation that addresses to identify how the porous media inside a TW can affect local flow phenomenon. This feature can be especially important to optimize porous media material, shape, and grain size to prevent dead zones and preferential flows, which can affect hydraulic retention time (HRT) and impact TW clogging and also affect biokinetics and treatment efficiency.

In this context, this study aims to propose a new approach using a CFD model coupled with a DEM model to provide a robust tool to support advanced modeling of TW combining fluid flow and solid-phase simulation. So, this article presents the CFD-DEM approach for studying the hydrodynamics of TWs and validates the results using experimental data from two case studies.

The proposed model was validated by comparing the simulation results with experimental data from two horizontal flow treatment wetlands (HFTW) case studies: Bandeiras (2009) and Gikas et al. (2017). It is important to highlight that calibration was not performed in the models because the aim was to verify if the features of the modeling approach enable it to provide results close to experimental data without the requirement of calibration. This premise was adopted to check if the modeling approach would be useful to model TW in future applications even if experimental data for calibration are not available.

The study by Bandeiras (2009) was conducted at the Laboratory of Environmental Sanitation (LSA) of the Department of Civil Engineering and Architecture (DECA) of the University of Beira Interior (UBI) in Portugal using a HFTW constructed in plexiglass, with a length of 2 m and a width of 0.80 m, filled with a 0.5 m high layer of expanded clay (diameter in the range of 4–8 mm).

According to Goniva et al. (2012), the coupling routine consists of defining the position and calculating the velocity of the particles by the DEM code, which is then sent to the CFD code. The CFD code determines the cell in the mesh to which the particle corresponds and then calculates the volumetric fraction of the particles and the average velocity of the particles. The fluid forces acting on each particle are then calculated, and the momentum between the particle and the fluid is defined by averaging all particles in the CFD cell. The fluid forces acting on the particle are sent to the DEM solver to be used in the next time step. Finally, the CFD solver calculates the fluid velocity considering the volumetric fraction of the particles in the cell.

The CFDEM^® coupling software enabled the coupling of CFD and DEM simulations. This software, developed by Goniva et al. (2012), is based on the open-source codes OpenFOAM (Open-Source Field Operation and Manipulation) and LIGGGHTS (LAMMPS (large-scale atomic/molecular massively parallel simulator) improved for general granular and granular heat transfer simulations), which are employed for fluid and granular particle simulations, respectively. The LIGGGHTS code is an improvement of LAMMPS.

The fluid phase was described using the Eulerian form of the Navier–Stokes equations and solved using an algorithm based on PISO (pressure-implicit with splitting of operators) algorithm in OpenFOAM, while the motion of the particles was described in the Lagrangian form of the equation of motion based on Newton's laws in LIGGGHTS code. The CFDEM^® coupling code was responsible for adding a volumetric fraction to the Navier–Stokes equations and exchanging information about the interaction between the fluid in the CFD component and particles in DEM.

Blender^® software was chosen for modeling the beds due to its open-source platform and ability to export 3D geometry in stereolithographic format (.STL), which is compatible with the software used in the model.

The simulations only included the active height of the beds. In addition, some simplifications were adopted for the inlet and outlet configurations. Indeed, the inlets and outlets were modeled as prismatic squared geometries instead of cylindrical pipes.

Among the solvers available in the CFDEM^® coupling software, the cfdemSolverPisoScalar was adopted, where the CFD component is based on the pisoFoam solver from OpenFOAM^®, which considers the flow of incompressible fluids in a transient regime and uses the PISO algorithm for resolution. Goniva et al. (2012) enabled the Lagrangian approach to the CFD solver and added a scalar transport equation.

Boundary conditions of the model were defined to indicate the flow inlet, outlet water surface, and wall conditions. For flow inlet, the massFlowInlet type was used, where the velocity is calculated based on volumetric flow, with the adopted value listed in Table 1. The boundary conditions adopted are summarized in Table 2.

A key aspect of CFD-DEM modeling is the mesh sensitivity analysis. In CFD-DEM models, the size of the CFD cell must be defined based on the particle diameter of the DEM. According to Peng et al. (2014), the size of the CFD mesh cell should be between two and four times larger than the particle diameter due to the nature of the algorithm in calculating porosity and drag force. However, excessively large cells can decrease the accuracy of CFD simulations. A mesh refinement was applied at the inlet and outlet zones to ensure better resolution in these areas.

Hence, mesh sensitivity monitoring for the different case studies was carried out with three different sizes based on the particle diameter. Pressure at a central point of the bed was used as the convergence criterion. The mesh was considered converged if the pressure difference at the point was less than 5%. Table 4 shows the details of the tested meshes for the different case studies.

DEM modeling is carried out by adjusting the input script, which is divided into initialization; definition of particle characteristics and insertion output options, and in case of coupling with CFD, specific commands must be inserted (Goniva et al. 2012).

Table 5 shows the main simulation settings. To enable information exchange between the simulation methods in CFDEM^® coupling, a command must be inserted in the input script of the DEM component authorizing the coupling. In Table 6, the main configurations for the coupling are presented.

The Bfull submodel was adopted due to its higher precision, as reported by Zhou et al. (2010). This submodel requires the use of the gradPForce force model, which calculates the force based on the pressure gradient. The force model is responsible for calculating the forces acting on each DEM particle and CFD cell. The VoidfractionModel was adopted due to the particles being smaller than the CFD mesh (Goniva et al. 2012).

Due to the lack of detailed information on the size distribution of the filling material in the case studies, the diameters of the aggregate type used in the experimental tests were 0.008 and 0.025 m for TW 1 and 2, respectively. Even if the use of only one homogeneous diameter is a simplification, it is expected that the fact of modeling the bed as a set of discrete elements in CFD-DEM will provide a model closer to reality than a continuous media used in the Darcy–Forchheimer approach.

For the expanded clay density in the simulation addressing TW 1, a value of 1,700 kg/m³ was adopted according to the manufacturer's specifications. For the density of the granite used in TW 2, the density of granite was considered at 2,600 kg/m³ (Olhoeft & Johnson 1989). As no studies were available regarding the friction of expanded clay, the values for natural rock were adopted. The rock–rock friction coefficient used was 0.808 (Dieterich 1972).

Table 7 summarizes the particle properties adopted in the study. Some values were adopted as default in the software.

In the simulation process, bed particles were inserted and settled under the influence of gravity. To save computational resources, the nve/sphere function was disabled after particle settling. This command is responsible for updating the position, velocity, and angular velocity of the particles, and disabling it avoids unnecessary recalculations since the particle positions have already been established. In addition, the time step was adjusted to the same value as the CFD component to ensure better integration between the simulation methods.

The coefficients were calculated based on the values of diameter and porosity provided by the authors listed in Table 1. The calculated D and F coefficients were 7,780,349.8 m⁻² and 2,640.60 m⁻¹, respectively, for the TW 1 and 3,386,718.8 m⁻² and 2,905.3 m⁻¹, respectively, for TW 2, using Equations (20) and (21). The pimpleFoam solver in the OpenFOAM^® software was used with the same parameters as in the CFD-DEM simulation.

In this study, concentration–time curves obtained from the CFD-DEM and Darcy–Forchheimer method were compared with experimental tracer data for TW 1 and TW 2, following the methodology described by Bandeiras (2009) and Gikas et al. (2017), respectively.

The concentration of the simulated tracer was measured using the probe function. This function collects the predefined data at each time step at a specific point in the simulation domain. In addition to the value of the passive scalar, the velocity and pressure values were also collected at each probe.

An error analysis of the simulated tracer HRT curve compared to the experimental tracer was performed based on the statistical metrics described by Fox (1981), Willmott (1982), Anthes (1983), Anthes et al. (1989), Pielke (2013), and Wilks (2011). Variance (s²), standard deviation (σ), mean error (ME), mean absolute error (MAE), mean squared error (MSE), root mean squared error (RMSE), RMSE_bias, concordance index (CI), and Pielke skill score index (D_PIELKE) was defined. The curves were compared up to the limit where the experimental measurements were made, i.e., 112 h for the Bandeiras (2009) case study and 197 h for Gikas et al. (2017).

The MSE (Equation (16)), which, despite being affected by singular points, is commonly used to measure the accuracy of numerical simulations.

The RMSE_bias (Equation (18)) removes a constant bias relative to the model's tendency from the RMSE.

D_PIELKE is an index proposed by Pielke (2013) that demonstrates the skill of a model based on the following criteria:

1. σs/σo close to 1;

2. RMSE < σo;

3. RMSE_bias < σo;

Tracer tests data from both case studies were used to validate the simulation models.

Monitoring the concentration–time curve of the tracer test allows the analysis of several aspects of the system's hydrodynamic behavior. In this study, the method of moments was used to analyze the tracer test response. In this method, it is necessary to determine the two moments of the concentration–time curve in the tracer test (Metcalf 2003):

Moment 1: the mean hydraulic detention time (M₁ or t_m), according to Equation (22), which indicates the average time that tracer particles remain in the system.

The contact time of the effluent with the bed elements is called the residence time, and in wastewater treatment systems, the longer the residence time, the longer the contact time of the effluent with the bed, and consequently, the better the pollutant removal efficiency (Danckwerts 1953; Levenspiel 2000).

RTD can also be analyzed based on the theoretical reactor models. There are two ideal flow models: the plug flow reactor (PFR) and the continuous stirred tank reactor (CSTR). In the PFR, fluid molecules cross the system in line and at the same velocity, so there is no mixing by axial dispersion. In complete mixing flow, the fluid composition is homogeneous due to the complete mixing of the liquid mass molecules (Danckwerts 1953; Levenspiel 2000).

The values indicated in Table 8 were used to classify the degree of the dispersion according to the δ value, based on Santamaria et al. (1999). The theoretical hydraulic detention time (τ) and hydraulic efficiency (λ) were also defined to diagnose the hydrodynamic conditions of the system.

In postprocessing, simulation results were compiled in the form of a color map for pressure and velocity, as well as velocity vectors. ParaView software was used for postprocessing. The pressure and velocity color map and velocity vectors can visually show the occurrence of deleterious processes, such as preferential pathways and dead zones, as well as demonstrate the phenomena occurring inside the bed (Fioreze & Mancuso 2019; Rengers et al. 2016).

This study is limited to the analysis of HFTW, but this approach can also be applied to other TW configurations. To simplify, the mesh simplifications were made in the geometry of inlet and outlet pipelines. So, the inlets and outlets were modeled as prismatic squared geometries instead of cylindrical pipes.

The software CFDEM^® coupling does not allow multiphase simulation, so the flow occurs in a saturated media. This assumption is common in TW models, as in the study by Ranieri et al. (2013). In addition, the software only considers the filling material particles to be spherical.

In Figure 3, concentration–time curves are presented, based on the monitoring of experimental and simulated tracers for TW 1 and TW 2, following the methodology described by Bandeiras (2009) and Gikas et al. (2017), respectively. These curves were compared with both the CFD-DEM proposed method and the classical Darcy–Forchheimer method.

In this study, concentration–time curves obtained from the CFD-DEM method were compared with the experimental data. The results indicate satisfactory agreement between the two approaches. For both case studies, tracer results obtained from the CFD-DEM model agree with the concentration peak obtained from experiments.

For TW 1 (Figure 3(a)), a difference of only 0.72 mg/L and 3 h in the peak of concentration was observed between experimental data and CFD-DEM simulations. The rest of the curve accurately matched the experimental data. For TW 2 (Figure 3(b)), the difference between simulation and experimental concentration peak was 1.58 mg/L and 0.34 days. It is important to note that in the study by Gikas et al. (2017), measurements were taken every 0.33 days, meaning the observed deviation differs in only one measurement step. This finding indicates that the flow velocity of the experiment and simulations using this approach agreed.

For TW 1, a disagreement occurred at the beginning of the curves. In addition, a long tail effect occurred, where small concentrations of the simulated tracer were detected for a long period after the saline tracer was not detected anymore by the measuring equipment. These findings indicate that the dispersion of the simulated tracer is lower than the experimental one in the first half of the curve but becomes adequate after the peak. Moreover, this elongation with small concentrations at the end of the curve may occur due to differences between the way the tracer was inserted in the experimental setup and in the simulation model.

However, in the experimental results for TW 1, a delay in the tracer's release was also identified, in which only 56% of the tracer was recovered. Bandeiras (2009) attributed this to the transport of tracer molecules in zones with little flow dynamics that may have subsequently evolved into dead zones, retaining the tracer inside these zones, as reported in the studies by Santamaria et al. (1999), Martinez & Wise (2003), and Albuquerque (2004). The long tail effect in the CFD-DEM simulation may corroborate this hypothesis, as the model may have underestimated the occurrence of dead zones in the flow, resulting in a slower release of the simulated tracer at the end of the curve.

Overall, the results obtained suggest that the CFD-DEM approach may be a viable alternative to simulate the transport of contaminants in packed bed systems. Although some differences were observed in relation to the experimental data, the simulated concentration values were very close to the measured values in both case studies. These findings indicated that the CFD-DEM model was capable of accurately capturing the main aspects of contaminant transport in TWs.

In addition, the analysis of concentration–time curves allowed the identification of some model limitations, such as the underestimation of tracer dispersion at the beginning of the curve and the occurrence of dead zones in the flow, which can lead to a long tail effect at the end of the curve. These limitations can be addressed in future studies through model refinement or the improvement of the experimental methodology, aiming to enhance the precision of the simulations and reduce the uncertainties in tracer analysis.

When comparing the experimental curves with Darcy–Forchheimer simulation results, several inconsistencies were observed between the simulation and experimental data. For TW 1, the simulated tracer curves were detected 1.08 days after the launch, with a peak concentration of 38.22 mg/L occurring after 2.38 days, a significantly higher and delayed value compared to the experimental data, where the peak occurred at 12.97 mg/L after 1.71 days of tracer launch. These inconsistencies demonstrate that the D and F coefficients failed to represent the behavior of the studied TWs. In TW 2, the detection of the synthetic tracer only occurred several hours after the experimental tracer had been completely extinguished, indicating that the simulation speed for this approach was much slower than the experimental speed, requiring coefficient calibration for a better curve fit. A long tail phenomenon also occurred in this approach.

Analyzing the tracer concentration curve distribution, it was observed that the CFD-DEM approach has a better capacity to indicate the behavior of the tracer inside the bed compared to the homogeneous porous media approach using Darcy–Forchheimer equations. It was noted that for both concentration–time curves using the classical approach, the shape is symmetric, similar to a normal distribution. In contrast, CFD-DEM results, especially TW 2, had an asymmetric shape, as observed in the experiment, indicating that the CFD-DEM approach is better at simulating flow than the Darcy–Forchheimer approach on this respect. This is due to the way the coupling between the codes physically inserts particles, better simulating flow in a granular media and being able to simulate small disturbances also detected in the experiment. On the other hand, the classical approach tries to simulate the porous media impact on the flow by attenuating the flow velocity, being incapable to represent local flow behaviors.

Table 9 presents statistical accuracy indexes for the two case studies, calculated by comparing the concentration–time curve values of the simulations with the tracer response in the experimental setups.

For TW 1, both CFD-DEM-simulated tracer curve and experimental tracer curve had similar deviations from the experimental results with a difference of only 1 h 30 min. In contrast, the Darcy–Forchheimer simulation had a deviation that was three times greater than the experimental results.

For TW 2, the Darcy–Forchheimer simulation came closer to the experimental results than CFD-DEM, but the difference between the latter model and the experimental results was only 0.79 days.

Concerning error indexes, the CFD-DEM simulations had presented substantially lower error indexes results than the Darcy–Forchheimer model. For TW 1, the MSE and ME were 30 and 27 times lower, respectively, whereas the MAE was six times lower. For TW 2, the MSE and ME were five and three times lower, respectively, whereas the MAE was two times lower for CFD-DEM.

Regarding CI results, the CFD-DEM model outperformed the Darcy–Forchheimer model in both case studies, with CI values of 0.92 and 0.76 for TW 1 and TW 2, respectively, compared to 0.29 and 0.00 for the Darcy–Forchheimer model. It is important to note that a CI of 1 indicates a perfect agreement between simulated and observed values.

The CFD-DEM model also demonstrated its proficiency in simulating TW according to Pielke's criteria, with D_PIELKE lower than the reference limit (D_PIELKE < 2) in both case studies. This is due to the small difference between the deviations of observed and simulated results, with both the RMSE and RMSE_bias being lower than the deviation of experimental results. However, the Darcy–Forchheimer model's suitability for TW simulation was not supported by D_PIELKE results. Indeed, in both case studies, the calculated values for D_PIELKE were higher than the aforementioned reference limit.

Table 10 shows the results of hydrodynamic parameters calculated using tracer studies and different simulation methods in the experiment. The table includes experimental data obtained from both Bandeiras (2009) and Gikas et al. (2017) and simulated tracer results for TW 1 and TW 2.

Comparing the values of t_m, in the simulation using the Darcy–Forchheimer approach, it was noted an interesting difference for TW 1, which is related to the discrepancy between the concentration–time curve of the simulation and the experimental curve. Indeed, this discrepancy led to t_m values for TW 1 and TW 2, which are 1.5 and 2.6 times higher than the experimental values, respectively.

In the CFD-DEM simulation, the difference between the t_m of the simulations and the experiment is 8 h 24 min for TW 1 and 2 h 10 min for TW 2. This difference can be attributed to the long tail effect, which tends to increase the t_m of the simulation, evidenced in TW1 whose curve had a longest tail in relation to TW2. Recalculating the parameters in TW1 without the tail extension, the hydraulic parameters tend to approach the experimental data. The recalculated t_m was 2.00 days for TW 1.

The long tail effect also affects the variance (σ²) of the CFD-DEM simulations, as it varies depending on the square of time, as presented in Equation (13). Removing the problematic part of the curve modifies the variance from 1.53 to 0.48 for TW 1 and from 3.21 to 2.71 for the TW 2 case study, whereas in the experiment, σ² was 0.63 and 2.05 for TW 1 and TW 2, respectively. The t_m value was higher than the τ value in both the experiment and in the simulations. Seeger et al. (2013) found values up to 50% higher than the theoretical value in their study and attributed it to tracer sorption in the support media or biofilm, but Bandeiras (2009) conducted an adsorption test in his study, ruling out the possibility of this phenomenon occurring in his experiment. This fact reinforces the possibility of stagnant zones retaining the tracer inside them.

The CFD-DEM simulation method can accurately determine hydraulic efficiency (λ) using Persson et al.’s (1999) parameter. In both case studies, the differences between the experimental and simulated parameters were minimal. For TW 1 experimental setup, λ value was 1.55, whereas the simulated value was 1.56. For TW2, the experimental λ was 0.97, whereas the simulated λ was 0.93. The similarity of results for λ from experimental and simulated values is directly linked to the fact that the CFD-DEM approach was able to precisely determine the peak time (t_p) of the tracer curve.

On the other hand, the λ values for Darcy–Forchheimer simulations were much higher than the experimental values: 2.36 and 3.15 for TW 1 and TW 2, respectively. This suggests that the homogeneous porous media approach is not as precise as the CFD-DEM approach in determining TW hydraulic efficiency.

Regarding the classification of the theoretical reactor flow model, TW 1 experimental data indicated a strong dispersion with N = 4, meaning that the flow slightly deviates from the plug flow. However, simulation results for the Darcy–Forchheimer model indicated considerably less dispersion, resulting in a N = 17, suggesting that this approach could not accurately classify the flow model. The CFD-DEM approach, in turn, classifies the flow as plug flow, with N = 6, which is only two tanks different from the experimental data. Indeed, the dispersion (d) from CFD-DEM results was slightly lower in the simulation than in the experimental study (Table 10).

In TW 2, experimental results also indicate that the flow approached to PFR, with N = 9, although there was still strong dispersion. The approach using the Darcy–Forchheimer equations also indicated flow close to plug flow, but with weak dispersion and N = 84, much higher than in the experiment (N = 9). In the CFD-DEM approach, the dispersion value is very close to the experimental value and the obtained N for the simulation is different in only 2 tank unit (N = 11) compared to the experimental result (N= 9).

In contrast, in the color maps for the simulation using the CFD-DEM approach (Figures 4(a) and 5(a)), the filler material particles are physically added to the simulation, which may cause small disturbances in the pressure along the bed. For TW 1, the pressure distribution pattern for CFD-DEM (Figure 4(a)) is similar to that of the homogeneous porous media (Figure 4(b)), with a gradual decrease along the length.

However, in the simulation for TW 2, there are significant differences in the pressure distribution pattern between the different simulation methods. In the simulation with the CFD-DEM coupling for TW 2 (Figure 5(a)), the pressure tends to decrease along the length, but with a heterogeneous distribution pattern, with several patches of higher and lower pressure along the profiles.

The difference in pressure distribution between the two CFD-DEM case studies can be attributed to differences in their characteristics, such as the inlet and outlet configuration, hydraulic loading rate, and physical characteristics of the bed, such as the aspect ratio, particularly the diameter of the packing particles. Indeed, in TW 1, both the distribution and collection of the effluent are spread along the width, making the pressure distribution pattern more homogeneous along the bed, while in TW 2, the collection of treated effluent is concentrated, causing pressure disturbances. The higher hydraulic loading rate in TW 2 can also be a factor contributing to the heterogeneous pressure behavior captured by the CFD-DEM model for this case study.

These findings indicate that the simulation method using the CFD-DEM approach allows for a better representation of the pressure distribution pattern compared to the simulation method of a homogeneous porous media using the Darcy–Forchheimer equations, given the inability of the classical method to describe differences in pressure distribution patterns.

For the simulation method using the Darcy–Forchheimer equations (Figures 6(b) and 7(b)), zones of higher velocity were observed near the inflow and outflow, and a constant and homogeneous pattern inside the bed. Stagnation zones were also formed at the corners. This behavior was also noticed by Rengers et al. (2016), who used the same homogenous porous media approach. Moreover, the velocity magnitude for the Darcy–Forchheimer method was lower than the CFD-DEM method, especially in TW 2 results, which is consistent with the results presented for the tracer concentration–time curves. For TW 1 using the classic approach (Figure 6(b)), only dead zones at the corners were identified. A slight velocity retardation was also identified at the bottom and sides, but it could not be considered a dead zone. Fioreze & Mancuso (2019), who adopted the Darcy equation for porous media flow using MODFLOW and MODPATH software, also identified a stagnation zone at the bottom of the TW in their simulations but did not make any considerations for the sides of the reactor. It is important to note that the CFD-DEM simulation method was capable to represent some hydraulic behaviors normally observed in wetlands that were not simulated by the Darcy–Forchheimer simulation method. For instance, CFD-DEM simulation results show that in both case studies, the flow velocity is higher at the top of the TW, indicating a possible preferential flow zone, and near the inlet and outlet zones. This behavior pattern was also observed in the experimental studies by Tanner & Sukias (1995) and Suliman et al. (2006), who identified an accumulation of solids near the bed surface, while De Paoli & Von Sperling (2013) noted solid accumulation near the inlet and outlet zones, indicating zones of higher flow rates.

CFD-DEM simulation results were also able to predict zones with low velocities on the walls and in the corners, pointing to possible dead zones. This behavior pattern was also found by Wang et al. (2014), where it was visually observed through fluorescent tracers for the configuration of inlet at the top and outlet at the bottom, similar to that used in TW 2, in which the flow passes via a short circuit at the top towards the outlet, forming dead zones at the bottom and in the corners of the wetland. Dead zones have a negative influence on the hydraulic efficiency of the wetlands. This reduction in effective volume is caused by several factors, including the relationship between length and width, filter material, and hydraulic application rate (García et al. 2005; Wang et al. 2014; Matos et al. 2018). Indeed, Fioreze & Mancuso (2019) recommend more studies related to flow inlet configurations in TWs, to prevent the formation of preferential flow zones and dead zones at the bottom of the system. The aforementioned study indicates the adoption of two distribution pipes at the inlet at different depths to avoid preferential paths and, especially, dead zones at the bottom of the bed. So, the capability of CFD-DEM to accurately identify the potential location of dead zones is useful information to support TW design and optimization.

The velocity vectors for the simulations using the Darcy–Forchheimer model (Figures 8(b) and 9(b)) move directly from the inlet to the outlet without anomalies in the flow direction, with a behavior similar to plug flow, as described by Fioreze & Mancuso (2019) in their model. In the velocity vectors using the coupling between CFD and DEM for TW 1 (Figure 8(a)), the behavior is similar to the simulations using the homogeneous porous media (Figure 8(b)).

For TW 2, small anomalies in the flow direction such as vortices and some vectors were detected in the CFD-DEM simulation (Figure 9(a)), especially near the inlet and outlet zones, probably due to the impact of particles on the flow in the zones with higher velocities.

On the other hand, the spatial distribution of the passive scalar for the CFD-DEM model indicated higher dispersion and areas where the tracer was retained, even after the majority had already exited the reactor. These areas support the zones of lower velocity detected in the velocity color map in Figures 6 and 7.

In fact, Bandeiras (2009) and Gikas et al. (2017), particularly for TW1, attributed the low tracer recovery obtained in their study to zones of low hydrodynamics that evolved into dead zones, retaining the saline tracer within. Therefore, the CFD-DEM model demonstrated its capability to accurately replicate the observed phenomenon, whereas the conventional model failed to capture such behavior.

Despite these simplifications, the simulation time for the CFD-DEM approach is much longer than for the simulation using the Darcy–Forchheimer model. The latter approach for simulating TW 1 case study takes only 15 min, compared to over 5 h required to the CFD-DEM to simulate 12 days of flow time. For the same flow time, the proposed approach for simulating the TW 2 case study takes over 20 h compared to about 15 min for the classical approach. So, the high computational demand of the CFD-DEM model limits the scalability of TWs analyses, requiring high-performance computers.

The present study aimed to experimentally validate a coupled CFD-DEM approach for modeling the hydrodynamics HFTW. In parallel, simulations using the classical homogenous porous media approach based on the Darcy–Forchheimer equations were also performed. Provided that no calibration was done, the comparison of simulation and experimental results showed that the CFD-DEM outperformed the Darcy–Forchheimer in terms of TW simulation based on statistical performance indexes and considering the similarity with experimental results. In addition, the CFD-DEM approach was better suited to represent details in hydrodynamics of the TW, including local flow anomalies and the internal behavior of the system. This capability of CFD-DEM is a relevant feature in simulations addressing biokinetic models or clogging assessments.

The application of CFD-DEM coupling for TW showed to be an interesting approach as it enables to perform more complex analyses of the hydrodynamics of TWs, becoming a valuable tool for future research and system optimization.

So, the proposed modeling approach was able to fill the objectives of the study opening the way for most advanced models using CFD-DEM, in terms of detailed assessment of local phenomenon as clogging simulation and preferential flow detailing. In future works, it is suggested to validate the model for other configurations of TWs and to incorporate additional processes, such as clogging and biomass growth. Although the CFD-DEM approach is more computationally expensive, its results show its potential to be applied in more in-depth studies of TW hydrodynamics.

The authors thank to Centro Federal de Educação Tecnológica de Minas Gerais (CEFET-MG) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) for the support.

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

The authors declare there is no conflict.