https://orcid.org/0000-0002-0840-6043
Abstract
Treatment wetlands are recognized as an effective technology for mitigating the impacts of urban runoff. However, there is no consensus on the design guidelines, and the effects of some design features, such as the underdrain system, remain unexplored. A simple analog model has been developed to mimic the underdrain network (when operating at saturation) and to evaluate the spatial heterogeneity of the flow entering it. The model has been applied to a treatment wetland in the Paris area and shows that the underdrain network was largely undersized, likely causing an uneven distribution of infiltrating flow along the length of the treatment wetland. It was also shown that this analog model can be used for optimization purposes and that it is important to use conservative values of the rugosity coefficient when designing an underdrain network.
HIGHLIGHTS
A brand-new analog model has been developed to mimic the hydraulic functioning of the underdrain.
Inflows into the underdrain decrease with the distance to the outlet.
Undersized underdrains can create uneven distribution of the percolating flow along the length of a treatment wetland (TW).
The analog model allows optimizing the design of underdrains for TWs with a saturated layer at the bottom.
INTRODUCTION
In recent decades, to attenuate the impacts of urbanization on the flow regimes and the water quality of receiving water bodies, various nature-based solutions have been developed, particularly the use of three types of vegetated systems: green roofs, bioretention (or biofiltration) systems and treatment wetlands (TWs). The latter two are similar in design. Bioretention systems are designed to store, filter and sometimes infiltrate urban runoff from an impervious surface at the source (Roy-Poirier et al. 2010) rather than being discharged through pipe systems to urban streams or other receiving water bodies. These systems are very effective in reducing urban runoff peak flows, stormwater volumes and pollutant loads (Liu et al. 2014). TWs, widely used for domestic wastewater treatment, have also become a popular solution for runoff treatment and retention (Shutes et al. 1999; Vymazal & Kröpfelová 2008; Choi et al. 2012; Schmitt et al. 2015; Adyel et al. 2016; Biswal et al. 2022) and combined sewer overflow treatment (Meyer et al. 2013; Pálfy et al. 2017; Tondera 2019; Rizzo et al. 2020). All studies converged on the same conclusion that TWs show promising performance for stormwater hydrology management (runoff volume and peak flow attenuation) and stormwater pollutant removal (total suspended solids, total nitrogen, total phosphorus and micropollutants) (Walaszek et al. 2018; Biswal et al. 2022).
Biswal et al. (2022) highlighted the need for guidelines for the design, installation, operation and maintenance of nature-based solutions to effectively remove environmental pollutants in the context of sustainable urban development. It is in this context of runoff water treatment by TW that the European Life ADSORB project takes place. One of its main objectives is to test the applicability of an innovative solution capable of efficiently reducing the loads of metallic and organic micropollutants, both in their particulate and dissolved phases. To achieve this goal, the project is based on the monitoring (operation and performance) of two TWs and their mechanistic modeling in order to improve the understanding of the removal mechanisms and to identify the most influential parameters for the water flow, transport and fate of micropollutants. The final objective of the Life ADSORB project is to determine the key points for the design of this type of TWs.
Shutes et al. (1999) attempted to provide some guidelines for the design, operation, maintenance and monitoring of a TW for the treatment of road runoff. They listed the factors that determine the design criteria for TWs, including traffic loads, road drainage area, land availability, cost and the size/extent and type of receiving water body. More recently, technical guidelines for the design, operation and management of constructed wetlands for stormwater treatment have been published in France (Molle et al. 2013) and Germany (DWA-A 178 2017). Biswal et al. (2022) indicated that the main operational and environmental parameters influencing the performance of TWs for stormwater remediation include hydraulic loading rate, media composition and depth, plant species diversity and seasonal variations. No study has addressed the design of the TW underdrain network, especially for long underdrains such as those in the Life ADSORB project. This study focuses on the underdrain network with the objective of evaluating the impact of its design on the flow distribution within the TWs. More specifically, this paper aims to answer the following two questions: (1) Is the flow rate uniformly distributed among the drain slots at constant load? (2) What is the effect of drain length and diameter on the spatial homogeneity of the percolated flows?
METHODS
Description of the studied site
The system studied in the Life ADSORB project consists in two 600 m2 vertical flow TWs located in the Bois de Boulogne park (Paris, France). Each was designed for the treatment of road runoff and, in particular, for the treatment of selected metallic and organic micropollutants. These two pilots are similar in terms of configuration (100 m long, 1 m deep) and operate with the same feeding period (alternating every month). The filter surfaces are flat. The only difference between the two is the composition of the filter layer, wherein one (TW1) consists only of sand (40 cm), while the other (TW2) consists of a specific layer of adsorbent material for the removal of micropollutants (Rainclean®, 20 cm) placed between two layers of sand (10 cm each). A 10-cm transition layer and a 50-cm drainage layer are present in both TWs. There is only one feeding point per TW, located at one end of each TW, and the treated water outlet is located at the opposite end. A throttle outlet located 30 cm above the bottom of each TW maintains a saturated layer and controls the outlet flow rate to not exceed 20 L/s when the TW is fully saturated. Since the configuration of the underdrain networks and the operation of the two TWs are similar, this study focuses on only one of the two TWs: TW1 and TW2.
Drainage system of the Life ADSORB TWs. The details of the underdrain network of TW1 are shown, while only the outline of the underdrain network of TW2 is shown in pink.
Drainage system of the Life ADSORB TWs. The details of the underdrain network of TW1 are shown, while only the outline of the underdrain network of TW2 is shown in pink.
The annual volume of water to be treated is estimated to be 200,000 m3/year. It is composed of runoff from the Paris ring road (21 ha), permanent parasitic clear water (279 m3/day) and combined sewer overflows (71 ha) in a minority way (a few percent of the annual volume). These waters are collected and stored in an underground structure with a capacity of 3,500 m3 and applied at controlled flow to the two TWs via a pumping station.
The maximum water level, called ‘high water level’, is set to 1.05 m above the bottom of TW1. A peculiarity of the studied system is that it receives water continuously (not only during storm events) due to the important parasitic water flow.
Evaluation of linear pressure head loss
is the linear hydraulic head loss in the main drain section 
, 
is the Darcy friction coefficient 
, L is the length of the main drain section 
, D is the diameter of the main drain section 
, V is the fluid velocity 
and g is the gravitational acceleration 
) is calculated differently depending on the value of the Reynolds number (
). If 
, it is estimated by Colebrook's equation:
and 
is the Reynolds number 
. Two values of the rugosity coefficient (
) have been chosen: 
m corresponding to a new PVC pipe and 
m corresponding to a PVC pipe with some deposit in it (Hager 2010).
is defined by
The linear hydraulic head loss in the main drain section is equal to 0.17 m for 
m and 0.21 m for 
m.
Therefore, the flow through the slots of the drain cannot be assumed to be uniform, and linear and singular hydraulic head losses in the drain must be considered.
Computation of slot singular pressure head loss coefficient
) was evaluated by computational fluid dynamics (CFD) using the COMSOL Multiphysics® software with a similar approach to that used by Forquet & Dufresne (2015). The computational domain is shown in Figure 2. It consists of a half-section of the drain (in blue in the figure) and a rectangular cuboid 0.06 m long (the distance between two slots plus the slot width), 0.5 m wide and 0.3 m high (the height of the saturation level) of the saturated porous media surrounding the drain (in gray in the figure).The slot, highlighted in blue in Figure 2, is the only connection between the two domains. Steady-state and incompressible Navier–Stokes equations were solved in the drain and slot domains, while Brinkman's equation with Forchheimer's correction was solved in the saturated porous media domain. The model was evaluated using nine combinations of boundary conditions (see Figure 2) ranging from 3,000 to 7,500 Pa for the pressure head above the drain (
), from 2,500 to 5,000 Pa for the pressure head downstream of the drain (
) and from 
to 
for the upstream flow rate (
). The range of 3,000–7,500 Pa for 
was chosen because it represents the amplitude of water pressure variation in the filter. The range of variation for 
was chosen based on the 
range to create sufficient gradients. Finally, the maximum flow rate at the surface is equal to the maximum flow rate that can be applied to the filter divided by first third of its surface area (because it is known that the distribution is uneven and occurs mostly in the first third of the filter). All combinations are given in Table 1. Subsequent mesh refinements were performed to achieve a reasonable mass balance. The resulting mesh consists of 943,363 elements. Finally, the singular head loss coefficient of the slot (
) is determined by fitting the relationship between the pressure head difference across the slot and the right-hand side of Equation (4). Boundary conditions and results of CFD simulations
| Pupstream (Pa) | Pdownstream (Pa) | Pressure head loss (Pa) | Upstream flow rate Qin (m3 h−1) | Slot flow rate (m3 h−1) |
|---|---|---|---|---|
| 3,000 | 2,500 | 500 | 0 | 0.19 |
| 3,000 | 2,500 | 500 | 36 | 0.18 |
| 3,000 | 2,500 | 500 | 72 | 0.16 |
| 7,500 | 5,000 | 2,500 | 0 | 0.44 |
| 7,500 | 5,000 | 2,500 | 36 | 0.43 |
| 7,500 | 5,000 | 2,500 | 72 | 0.41 |
| 7,500 | 2,500 | 5,000 | 0 | 0.62 |
| 7,500 | 2,500 | 5,000 | 36 | 0.62 |
| 7,500 | 2,500 | 5,000 | 72 | 0.60 |
| Pupstream (Pa) | Pdownstream (Pa) | Pressure head loss (Pa) | Upstream flow rate Qin (m3 h−1) | Slot flow rate (m3 h−1) |
|---|---|---|---|---|
| 3,000 | 2,500 | 500 | 0 | 0.19 |
| 3,000 | 2,500 | 500 | 36 | 0.18 |
| 3,000 | 2,500 | 500 | 72 | 0.16 |
| 7,500 | 5,000 | 2,500 | 0 | 0.44 |
| 7,500 | 5,000 | 2,500 | 36 | 0.43 |
| 7,500 | 5,000 | 2,500 | 72 | 0.41 |
| 7,500 | 2,500 | 5,000 | 0 | 0.62 |
| 7,500 | 2,500 | 5,000 | 36 | 0.62 |
| 7,500 | 2,500 | 5,000 | 72 | 0.60 |
CFD geometry used to estimate the singular head loss of the slot. Boundary conditions are also shown (red areas for pressure boundary conditions in the porous media (
) and downstream (
) of the drain, and arrow for flow boundary conditions (
)).
CFD geometry used to estimate the singular head loss of the slot. Boundary conditions are also shown (red areas for pressure boundary conditions in the porous media () and downstream () of the drain, and arrow for flow boundary conditions ()).
The analog model
An analog model of the drainage has been developed that considers the effect of the drainage network geometry (pipe diameter and antennas) and linear and singular head losses on the distribution of flow through each slot of the drainage network.
in the underdrain at slot j is set to a constant value 
. Then the algorithm calculates the inflow through each slot from upstream to downstream of the underdrain network. These values are used to calculate the flow along the entire underdrain 
to the outlet. Knowing the flow rate 
, the linear head losses between the slots are then estimated. Finally, 
are updated based on the linear head losses starting downstream. The loop continues until the differences between the slot level pressures 
obtained in the previous iteration and those obtained in the current iteration (Delta_P) are less than a given threshold (tolerance). 
Analog model flow chart.
The calculations performed at each step of the algorithm are described in more details hereafter:
- The flow rate passing through each slot is calculated using Equation (4):where(4)
is the singular hydraulic head loss 
, K is the head loss coefficient of the slot 
, 
is the flow rate passing through the slot 
and 
is the cross-sectional area of the slot 
. The estimation of the head loss coefficient (
) is discussed in detail further in the Methods section. The flow rates passing through the pipes between each slot are calculated iteratively from upstream to downstream by summing the flow rates passing through the slots located upstream.
The linear hydraulic head losses associated with the flow rates calculated in the previous step are calculated (one value for each section of a pipe located between two slots).
- The pressure head at the outlet of the underdrain network (
) is calculated by estimating the singular head losses associated to the cumulative flow rate passing through the throttle outflow,where(5)
is the flow rate exiting the drain (calculated by summing up the incoming flow rate passing through each slot) (
), 
is the singular head loss coefficient (equal to 0.62 for a circular opening) 
and 
is the height of the throttle opening measured from the bottom of the TW 
. The pressure heads at the level of each slot are obtained iteratively starting from downstream to upstream by adding to the pressure head at the slot downstream the linear head losses associated with the flow through the section of pipe located between the current slot and the slot downstream.
- To ensure model convergence, a regularization step has been introduced. The pressure heads at each slot (obtained from the previous step) at iteration i are adjusted following the equationwhere(6)
is a regularization parameter comprised within 0 and 1 (0.99 for this study). - Finally, the pressure heads at each slot at iteration i are compared to pressure heads at iteration
and the loop can be exit whenwhere(7)
is the tolerance, set to 
in this study.
Sizing optimization algorithm
) for the slot, the analog model was used to optimize the diameter of the drain in order to minimize the difference in flow rate between the slots. It is assumed that the diameter of the two antennas must be identical. Therefore, two diameters have to be optimized: the diameter of the main drain and the diameter of the antennas. For this purpose, an objective function (
) has been defined that measures the difference in flow rate between the first and the last slot of each branch,
being the flow rate passing through the first slot of a branch (
) and 
being the flow rate passing through the last slot of a branch (
).The objective function is evaluated independently for the main drain and the two antennas. The results are stored in a triplet. Starting from the existing diameter, the diameters of the main drain and the antennas are iteratively increased if they are above a fixed threshold (
in this study) until all values of the triplet are below the threshold. The process of increasing the pipe diameter must satisfy two conditions: (i) the diameter must be a commercial diameter and (ii) the main drain diameter should be greater than or equal to the antenna diameter.
RESULTS AND DISCUSSION
Slot singular pressure head loss coefficient
183 
. The resulting value is in the upper range of values found in the literature but is expected considering the small size of the slot and the fact that it acts as a diffuser leading to significant energy losses (Hager 2010). 
Determination of the singular head loss of the slot by fitting results of numerical experiments.
Determination of the singular head loss of the slot by fitting results of numerical experiments.
Estimation of flow distribution along the underdrain network using the analog model
) ranging from 0.4 to 1 m. The resulting distribution of flow rates across the slots is shown in Figure 5. Most of the linear head losses occur in the main drainage pipe, resulting in a flow rate divided by nearly 3 between its outlet and inlet. High linear head losses require higher upstream pressure to allow the water to flow, and therefore reduced pressure head differences between the surface and the pressure head in the pipe for the upstream slots. It is also observed that the outflow rates in the antennas vary significantly. Overall, the results of the analog model suggest that the underdrain network was largely undersized. The consequence of the heterogeneous flow rates through the slots is the creation of a dead zone upstream and, consequently, the development of a preferential path that shortcuts part of the porous media. 
Distribution of simulated flow rates across each slot of the underdrain network as a function of its distance to the underdrain network outlet for different water heights in TW1 (rugosity 
).
Distribution of simulated flow rates across each slot of the underdrain network as a function of its distance to the underdrain network outlet for different water heights in TW1 (rugosity ).
Optimized sizing of the underdrain network
In order to reduce the linear head losses, pipe diameters were optimized using the algorithm described in “Sizing optimisation algorithm” paragraph in the “Methods” section. Calculations were performed for pipes with smooth and rough rugosity (
and 
, respectively), and the results are presented in Table 2. Optimized values are significantly larger than the current design values, indicating that the underdrain network was not designed to distribute the inflow rate uniformly along its length. The same optimized diameters were obtained under the smooth and coarse assumptions, and the optimized diameter for the main section is larger than that of the antennas.
Drainage pipe diameters in mm: current values versus optimized ones for smooth and rough pipes
| Pipe | Current diameter (mm) | Optimized diameter (mm) (k = 1.5 × 10−6 m) | Optimized diameter (mm) (k = 1 × 10−4 m) |
|---|---|---|---|
| Main | 160 | 315 | 315 |
| Antennas | 160 | 250 | 250 |
| Pipe | Current diameter (mm) | Optimized diameter (mm) (k = 1.5 × 10−6 m) | Optimized diameter (mm) (k = 1 × 10−4 m) |
|---|---|---|---|
| Main | 160 | 315 | 315 |
| Antennas | 160 | 250 | 250 |
. It is plotted on the same y-axis as Figure 5 to illustrate how reducing the linear head losses allows for a more homogeneous distribution of flow through the drainage pipe slots. The results are not significantly affected by the height of the saturation zone in the TW. 
Distribution of simulated flow rates across each slot of the underdrain network as a function of its distance to the underdrain network outlet for different water heights in TW1 (rugosity 
).
Distribution of simulated flow rates across each slot of the underdrain network as a function of its distance to the underdrain network outlet for different water heights in TW1 (rugosity ).
CONCLUSIONS
After identifying a problem with the flow distribution in the underdrain network using a simple pressure head loss calculation, a simple algorithm was developed to estimate the flow through the slots of the underdrain network under constant loads. The resulting model shows the following under the current design:
The flow is not uniformly distributed along the slots.
The flow is stronger for the slots located near the outlet and the difference increases with the water height in the TW. For a water height of 1 m, close to the ‘high water level’, the flow at the last slot reaches 0.6 L/s, while for a water height of 0.4 m, no flow passes through the first slot at the beginning of the longest antenna.
The current design is likely to create preferential flow paths and dead zones within the filter media, thereby reducing the active volume of the filter media and the water residence time, thereby compromising the treatment capabilities.
The developed model can be used not only to quantify the spatial heterogeneity of flow in an underdrain network but also to identify where most of the head losses occur.
The main drainage section was identified as the bottleneck of the current design, as its diameter is largely undersized.
The model was successfully used to optimize the sizing without changing the network layout to reduce the spatial heterogeneity of the flow.
It also demonstrated the importance of anticipating the effects of pipe fouling by using a rugosity coefficient representative of a mature system.
This study demonstrates the importance of careful design of the underdrain network, which is often overlooked in the field of TWs. A numerical model has been developed for this purpose, but it is currently limited to TWs with a permanent saturated layer at the bottom (where the underdrain network is located). The model presented in this paper is primarily for design purposes. For applications where a dynamic approach is required, it can be coupled with a mechanistic numerical model of the TW. In particular, this will allow the user to consider cases where the water level is not the same along the entire length of the system. Finally, the optimization step could be further improved not only by adjusting the pipe size but also by identifying the best layout for the underdrain network.
ACKNOWLEDGEMENTS
This project has received funding from the European Life ADSORB project (LIFE17 ENV/FR/000398).
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
CONFLICT OF INTEREST
The authors declare there is no conflict.
%20cropped.png?versionId=5945)