Hydraulic simulation of perforated pipe systems feeding vertical flow constructed wetlands

Constructed wetlands are a combination of soil filters (porous media such as sand or gravel) and marsh plants such as the common reed, Phragmites australis. The subsurface flow of wastewater is either vertical or horizontal. Whereas vertical flow constructed wetlands are fed intermittently and therefore show varying saturation levels in the soil filter, horizontal flow constructed wetlands are fed continuously and work under fully saturated conditions (Giraldi & Iannelli 2009). In both systems, pollutants are removed via biological, physical and chemical processes during the underground passage (Aiello et al. 2016).

Vertical subsurface flow constructed wetlands (VFCWs) have become an accepted alternative means for the treatment of wastewater of different origins such as households, dairy production, industries, combined sewer overflows, storm water, etc. Most research has focused on treatment performance (Al-Rubaei et al. 2016), filter layout (Zhao et al. 2004), solids and biofilm accumulation (clogging) (Zhao et al. 2009), and hydraulics in the filter (Langergraber et al. 2009; Meyer et al. 2015). An aspect that is mostly overlooked is the distribution of wastewater on the surface of VFCWs. A uniform distribution is essential for the treatment of wastewater, as highly loaded parts of the filter are prone to clogging and do not remove pollutants reliably; filter space with low hydraulic load, on the other hand, does not contribute to the cleaning process and is therefore inefficient and economically unsound (Persson et al. 1999).

Wastewater distribution on the surface of VFCWs is accomplished via porous pipe systems made of plastic. Few rules are found in the technical literature as to the design of such pipe systems. The rules of technology in Germany (DWA 2017) require that openings be circular, with a diameter of no less than 8 mm and a maximum area that is supplied per opening of 5 m². The feeding flow rate in the inlet pipe has to rise up to 6 L/(min m²) or more. As rules are not very detailed, most practitioners rely on their experience when designing pipe systems for the feeding of VFCWs. For pipe systems in operation, an evaluation of the spatial distribution of the discharge is virtually never undertaken.

The aim of this study is to investigate whether hydraulic modeling is a feasible and reliable method for assessing the uniformity of discharge from perforated pipe systems. First, the requirements and conditions for the simulation of perforated pipe systems are investigated and the model basics are explained. Then the whole process of model build-up, calibration, application and analysis is presented and discussed. For the hydraulic computations, the software EPANET is employed. Pressure measurements in the pipe systems are used to assess the performance of the VFCWs and provide data for calibrating the model. Calibration is done via a genetic optimization algorithm called NSGA-II. Through calibrated model results, the aging (clogging) stage of a given wetland and the discharge distribution of its pipe system are assessed.

Scientifically, the modeling of perforated pipes is related to the simulation of irrigation pipes (agriculture) and of manifolds for fluid distribution or heat transfer in a number of technical applications, e.g. fuel cells, heat exchangers, spargers or solar collectors (Chen & Sparrow 2009; Wang 2011). Wang (2011) names three approaches for the study of flow distribution in manifolds: discrete (network) models, analytical (continuous) models, and computational fluid dynamics (CFD).

For its relative simplicity and therefore easy handling, this study uses a discrete (network) model to describe flow patterns in and discharge from perforated pipe systems. According to Smith et al. (1986), four equations are sufficient for this purpose.

• (1)
  Mass continuity for any given control volume:

  

  (1)

  where Q is the flow rate.
• (2)
  Energy conservation:

  

  (2)

  where z is the elevation head; p represents the pressure; ρ is density; g is gravity acceleration; v is the mean velocity of the fluid; h_v is total head loss.
• (3)
  Pipe friction:

  

  (3)

  where head loss (h_v) is computed as the sum of friction losses (h_r) and local losses (h_v,l). Explaining the details of calculating friction losses is beyond the scope of this paper. Details can be found in the literature on hydrodynamics (e.g. Crowe et al. 2005; Krause 2005). Local losses will be discussed further in the Local losses section.
• (4)

  Outflow characteristics of openings: see Outflow characteristics section. Pipe networks that feed VFCWs are often non-trivial in their layout. Commonly, a distributor pipe is attached to an inlet pipe. From the distributor, branches with small holes in them stretch over the soil filter (see Figure 2). The links between distributor and branches are realized via tee junctions (dividing junctions). For diameter reduction between inlet and distributor pipe as well as between distributor and branch pipes, contractions are installed. To take these geometric features into account, we add a fifth requirement to the simulation.

• (5)

  Local losses in pipe fittings such as bends, contractions, T-shaped fittings (tees), etc.: see Local losses section.

It is known that with the aging or clogging of a wetland, hydraulic conductivity is reduced in the soil filter (Zhao et al. 2009). For the model we assume that reduced hydraulic conductivity affects the outflow characteristics of the pipe system: lower conductivity results in a kind of flow resistance or back pressure that leads to reduced outflow from the openings. This reduction in discharge can be represented by the discharge coefficient μ. In that way, determining μ by calibrating a model can provide information about the state of clogging in a given wetland.

The study was undertaken on a privately owned VFCW of the size of 25 population equivalents. In the wetland, the domestic sewage of a small farm (inhabitants, workers, visitors) is treated. The wetland has a surface of 114 m² and is fed by two porous pipe systems (Figure 2), which are covered with a 5 cm layer of gravel. The systems are fed in parallel from the same primary treatment tank. The feeding is volume-based rather than time-based. Treatment performance of the wetland is good (mean value chemical oxygen demand (COD) effluent: 35.2 mg/L, efficiency: 82%).

For the measurement period of about 12 weeks, the entrance to one of the pipe systems was blocked and only the second pipe system (feeding a surface of 60 m²) was supplied with water (Figure 2). This was necessary because feeding both systems in parallel did not produce pressure data that could be properly analyzed. A feed line of 19.85 m length, a distributor pipe of 6 m length, and 19 branch pipes of around 10 m length make up the perforated pipe system. Diameters are 100, 46, and 28 mm, respectively. Every 25 cm, the branches are drilled with circular openings of 3 mm diameter. The system is fed by gravity from a siphon with an elevation of around 1.2 to 1.33 m (depending on water level) over the invert of the branch pipes; no pump is employed. Due to imprecise work during construction, the filter surface is uneven, with as much as 30 cm difference in elevation between far ends. These differences in elevation were implemented in the model as well (pipe openings).

Six pressure probes of type PR36XW by the firm Keller were installed in opposite ends of three of the branch pipes and another probe at the end of the inflow pipe, right before the first tee junction (Figure 2). The data were logged second by second onto loggers of the firm Hobo. From the measurements, it can be concluded that each feeding interval consists of three phases (Figure 3). First, when water enters the pipe system, a pressure rise is seen. Second, a stationary phase is established when the pipes are fully filled with water. Third, the inflow from the tank is stopped; while the remaining water leaves the pipe system, a pressure drop is registered. The arithmetic mean of the pressure values during the stationary phase is taken as the representative pressure value of that interval. In this paper, results for measurements from 10 October through 1 November 2016 are presented. To minimize exterior disturbances, intervals on rainy days and with less than 6 hours' percolation time between feedings were excluded from the analysis. The mean values of the remaining values were then taken as reference values for the calibration process.

The pipe network of the wetland was created using the software environment R (R Core Team 2017). For the hydraulic computations, the program EPANET (https://www.epa.gov/WATER-RESEARCH/EPANET; Rossman 2000) was employed. EPANET is a public domain software provided by the United States Environmental Protection Agency and is used world-wide for the modeling of pressurized drinking water networks.

For each pipe, local loss coefficients can be specified in EPANET. For bends and contractions, values from the literature (Idelchik 1996) were taken. As loss coefficients in tee junctions depend on flow, an iteration sequence was programmed in R. The model is fed with start values for the loss coefficients. From the model results, loss coefficients are calculated according to Equation (6) and Table 1 and fed to the model. Again, the results are taken to calculate the loss coefficients. The iteration stops when the deviation of the sum of all the loss coefficients from one iteration step to the next falls below a given threshold (0.1).

Smaller values represent higher uniformity of emission rates. Values under 5% are classified as good, between 5% and 10% as medium, greater than 10% as poor (Hezarjaribi et al. 2008).

Before optimization, a comparison of the different approaches for calculating loss coefficients in tee junctions (Table 1) was undertaken. Five versions of the hydraulic model were run, four with iterative calculation of the coefficients, one with constant values taken from Rossman (2000).

Calibration is done by an automated optimization algorithm. Thereby, tedious ‘trial-and-error’ calibration can be avoided, the calibration process can be accelerated, and goodness of fit can be improved and modeler bias reduced (Dent et al. 2004). Having two objective functions (RMSE and Dev_Q) calls for an algorithm that can do multicriterion optimization. In this study, the NSGA-II (nondominated sorting genetic algorithm II) is employed (Deb et al. 2002). The algorithm successively takes samples from the search space (called populations), evaluates them and evolves them into a new population (the next generation). Optimization stops after a given number of generations of a given size. The NSGA-II does not find a single solution but a set of Pareto-optimal solutions. When comparing any two Pareto-optimal solutions, none of them can be considered ‘better’, as one will yield a better result in terms of one performance index (e.g. RMSE) and the other will provide a better fit in terms of the other performance index (e.g. Dev_Q). In the end, the modeler has to select from the Pareto set that solution which is most suitable for their optimization problem. In this study, Pareto-optimal solutions with a Dev_Q < 3% (an achievable value according to preliminary model runs) are filtered, and of the remaining solutions, the one with the lowest RMSE is selected.

Model parameters μ (discharge coefficient) and ζ_inlet (the combined loss coefficient of the inlet pipe with pipe entrance in the tank and fittings of unknown number and geometry) are selected as calibration parameters. Calibration is done in two stages (Figure 4), where the first stage serves as a test for parameter limits and the computation method of loss coefficients. Computation after Gilman (1955) was not taken into account, as coefficient values are close to those by Idelchik (1996), so that results can be regarded as equivalent. Population size and generations for the NSGA-II are set in such a way that computational cost remains low while optimization shows good convergence. These values as well as initial parameter values were chosen after preliminary calibration runs. For stage two, one method for computation of loss coefficients in tees is chosen, and parameter limits are set in such a way that the space is only slightly larger than the space that the Pareto set from stage one occupies. By setting smaller parameter limits, it is assumed that a better solution (in terms of better fit) is accomplished.

The arithmetic means of the representative pressure head values in the pipe system are between 0.222 and 0.475 m for the seven probes. During a second measuring campaign, comparable values were achieved (maximum difference of 22 mm), supporting the plausibility of the measurements. From the measurements in the tank, a mean feeding volume of 486.9 L per interval is calculated, while the outflow during the stationary phase comes to 2.811 L/s (2.861 L/s in the second campaign). Feeding intervals are registered every 4.7 h, resulting in a daily mean hydraulic load on the wetland of 41.2 L/(m² d). Although percolation time is less than the 6 h prescribed in the technical rules in Germany (DWA 2017), the hydraulic load is considerably less than the allowed maximum of 80 L/(m² d). Therefore, the wetland is supposedly not overloaded from the hydraulic point of view, and treatment performance should be ensured. The feeding flow rate is only 2.80 L/(min m²), less than half the flow rate of 6 L/(min m²) prescribed in the technical rules.

Local loss coefficients that are calculated according to different approaches are compared in Figure 5. For coefficients in the lateral flow direction, higher values (up to 81.13) result for tee junctions closer to the inlet of the distributing pipe (tee #1 and following). This means that water that branches off there receives comparatively higher losses than water that branches off later. The values are lowest for Winter (1955) and highest for Gilman (1955) and Idelchik (1996). Straight direction flow loss coefficients are 0 at the beginning of the distributor for Idelchik (1996), Gilman (1955) and Glück (1988) and rise up to no more than 0.36 for Glück (1988), with other values below that. Loss coefficients after Winter (1955) are 0.1 and higher, whereas loss coefficients according to Idelchik (1996) stay 0 in this case for all straight direction tee ends.

The difference the loss coefficients make in the distribution of discharge from the pipe openings can be seen in Figure 6. The best uniformity is achieved with Winter (1955) (CV_Q = 12.5%). Boxplots after Gilman (1955), Glück (1988) and Idelchik (1996) (CV_Q = 24.9, 21.5 and 24.9%, respectively) are a little taller, but still small when compared with that which shows the discharge with loss coefficients taken from Rossman (2000) (fixed values of 0.6 for straight direction flow, 1.8 for lateral direction flow). Here, the CV_Q goes up to 59.8%. In some of the openings, negative discharges occur, a very unrealistic scenario.

From these results it is concluded that selecting loss coefficients should be done with great care, as they strongly influence model outcome. Implementing constant values is discouraged. The comparison does not give a hint as to which is the best method for calculating loss coefficients. Gilman (1955), Idelchik (1996) and Glück (1988) produce results that are closer to each other than when compared with values from Winter (1955).

One drawback in using EPANET for hydraulic modeling of small pipe systems is that it does not accept negative values for loss coefficients; therefore, negative values were replaced by 0 when fed to the model. Supposedly, this has an effect on model results such as flow patterns and discharge from the pipe system. Negative loss coefficients occur in the straight direction branch of tees when pressure rises after flow branching, as is commented on in a number of publications (e.g. Wang 2011; Idelchik 1996). The reason is that the fluid that branches off has a lower energy level (lower velocity) than the fluid in the pipe center (higher velocity). When lower-energy fluid leaves the main pipe, the average energy level rises, which can be expressed by negative loss coefficients. Coefficient computation after Idelchik (1996) results in negative coefficients for v_s/v > 0.5, but a value of 0 has to be implemented in this study.

The calibration results chosen out of the Pareto sets in stage one are summarized in Table 2. Calibrated parameter values for Glück (1988) and Winter (1955) are close to each other, but Glück (1988) results in a better goodness of fit. Parameter values for Idelchik (1996) show a slightly higher μ and smaller ζ_inlet and result in the lowest RMSE. Therefore, in the second stage, calibration loss coefficients in tees were calculated after Idelchik (1996).

Stage two resulted in a set of eight Pareto-optimal solutions (Table 3), all of which fall below the chosen limit for Dev_Q of 3%. In terms of goodness of fit, the results are acceptable. Calibration did not provide one ‘perfect’ set of values. According to the rule established earlier (parameter set with Dev_Q ≤ 3% and lowest RMSE), the solution μ = 0.215 and ζ_inlet = 94.5 is selected.

A model can never be a perfect representation of reality, so discrepancies between measured and calculated values can be explained by inadequacies in model structure as well as model parameter values. Also, measuring errors influence the model fit.

From the calibrated value of μ = 0.215 we can conclude that the perforated pipe system and/or the soil filter are significantly aged, as the value is much lower than the standard value between 0.6 and 0.7 and even below the range of 0.4 to 0.9 given by Smith et al. (1986). As the pipe system is regularly maintained and cleaned (high-pressure flushing once a year with 80 to 100 bar) and is known to have no defect, the effect will most probably be caused not in the pipe system but in the filter, where after eight years of operation the soil pores must be partly clogged by suspended solids and biofilm. Year-by-year assessment of effluent concentrations shows that treatment performance, however, seems not to be affected negatively by the clogging. Some sources (e.g. Al-Rubaei et al. 2016) assume that a certain stage of clogging is even beneficial for the removal of pollutants, as the biofilm that causes the clogging holds microorganisms that metabolize organic substances.

The wastewater distribution according to the hydraulic simulation is seen in Figure 7. Discharge uniformity is judged as poor, with CV_Q as high as 11.7%. The maximum discharge per opening (Sector A) is about 60% higher than the minimum (Sector B). Assuming pollutants to be evenly dispersed in the wastewater, sector A receives a higher load of COD and NH₄-N than sector B. Possibly, clogging is stronger in sector A than in other parts. If clogging reaches a critical state, visible results such as ponding will occur first in that area. However, the mean daily load in the heavily loaded area is only 7.9 g COD/(m² d), way below the maximum load of 20 g COD/(m² d) prescribed in the technical rules in Germany (DWA 2017).

We assumed that μ takes a constant value all over the surface. This assumption might not be valid: if hydraulic loads are non-uniform, non-uniform clogging in the filter can occur. This will result in higher μ values for areas with less loading and vice versa, an effect that was not taken into account in this model. On the other hand, areas with higher clogging and smaller μ will receive less hydraulic loading, which might help to reduce clogging in those areas and lead to more uniform μ values on the filter surface in the long run.

Model optimization strongly depends on the sophistication of the model in structure and parameters (Dent et al. 2004; Savic et al. 2009), so geometry and parameters were implemented with great care. However, any model can only be an approximation to reality. As for parameter μ, the implementation as a constant factor for all openings is questionable due to drilling burrs, slime clinging to holes, or plant growth that might obstruct outflow differently in different places. Secondly, all these phenomena, not just clogging, could influence the value of μ. Since the pipes are regularly maintained with high-pressure flushing, we assume the openings are not obstructed, so the simplification of a constant factor is acceptable. And, although we cannot totally eliminate other causes than clogging for the resulting value of μ, given the good maintenance of the pipe system, we attribute the biggest part of the effects we found to clogging. Moreover, it would be unfeasible to calibrate μ for every single one of the 741 openings in the model, so some level of simplification is unavoidable.

Validation and verification of model assumptions and results would be desirable. Measuring water distribution in the soil filter could be employed to confirm the modeled wastewater distribution on the surface. Giraldi & Iannelli (2009) present a methodology for measuring water content profiles via capacitance probes. So far, little is known about the link between discharge from a pipe system and the hydraulic conductivity in the filter or the state of clogging. A study on a young wetland spanning several years, employing repeated pressure measurements and calibration of μ to detect aging, and possibly accompanied by clogging assessment (methodology overview in Nivala et al. 2012), can help to detect such a link. Laboratory scale tests on columns with accelerated clogging would also be helpful. Also, the methodology presented in this study should be employed for other wetlands and/or the same wetland in the future to confirm if preconditions and assumptions hold and if results are robust.

Essentially, modeling pipe systems can be thought of as a prerequisite for simulating species fate in the soil filter: it is hydraulics that defines the reaction room, and differences in discharge on the filter surface will lead to different reaction rates and results in different parts of the soil filter.

For a VFCW in operation, loadings can be reviewed via pipe system models. One aspect is the uniformity of distribution of wastewater on the filter surface. Another aspect is to check whether the volume given onto a surface during one feeding is sufficient for uniform distribution. If needed, restoration strategies for clogged soil filters can be construed from that. Modeling perforated pipe systems can be a helpful tool for engineers when designing such networks. Modeling allows quick comparison of different network layouts, and the cost for building such models can be cut enormously by automatically building network representations in computer code. By implementing different pumps, performance and costs can be evaluated. Aging of a pipe network and its consequences for flow patterns can be estimated by adjusting the value for μ in the model, allowing an estimate of network performance in the future.

The model presented here only accounts for the stationary phase of a feeding interval. No approaches for modeling discharge during the phases of pressure rise and pressure drop have been formulated. It can be assumed that, depending on the size of the pipe system, a certain amount of water is discharged from the openings during these phases as well. For a more detailed assessment of the wastewater distribution, future models should take that into account.

• EPANET, software that was developed for simulating drinking water networks, can be employed for the simulation of small perforated pipe networks and the evaluation of their performance (i.e. the distribution of wastewater on the surface of VFCWs).

• For hydraulic modeling of small pipe systems, local loss coefficients ζ for dividing junctions are to be implemented with great care, as they strongly influence model outcomes. Calculating loss coefficients according to flow velocities in the pipe system is encouraged.

• Obstructions to discharge from pipe openings, such as clogging in the soil filter, drilling burrs, and plant growth, feed back on the pipe system model via the discharge coefficient μ. Calibrating μ on a wetland with a well-maintained pipe system allows an assessment of the clogging stage of the filter, as other effects are assumed to be negligible in that case.

• For calibration, a multicriterion algorithm is useful as it can handle multiple objectives at once. From the resulting set of solutions, the modeler has to select the most appropriate solution for their purpose.

• Through hydraulic simulation, the distribution of wastewater on a filter through a porous pipe system on the surface can be assessed. Modeling can also be employed in the designing stage of a network to evaluate distribution performance, pumping regimes, etc.

• Pressure measurements in the pipe system are advisable in order to support hydraulic simulation of such a system and to characterize the wetland that it feeds.

• This study presents a simple procedure to create significant findings on wastewater distribution on the surface of a VFCW. Sole preconditions are geometric information on the feeding pipe system, pressure measurements in the pipes, and information on the pump or tank feeding the pipes.

The study that is presented in this paper was funded by the German Ministry of Education and Research via the project management agency Karlsruhe, Department Water Technology, support code 02WQ1302. Also we acknowledge the support of farm Hof Mahlitzsch. The said sponsors had no involvement in the study design, the collection and interpretation of data, or the writing of the report.