# Sizing efficient underdrains for treatment wetlands

Ania Morvannou, Matthieu Dufresne, Marie-Christine Gromaire, Stéphane

Troesch, Nicolas Forquet

<figure><img src="image_1.png" /></figure>

## To cite this version:

Ania Morvannou, Matthieu Dufresne, Marie-Christine Gromaire, Stéphane Troesch, Nicolas Forquet. Sizing efficient underdrains for treatment wetlands. Water Science and Technology, 2023, 89 (2), pp.241 - 251. 10.2166/wst.2023.417. hal-04564889

HAL Id: hal-04564889

https://hal.inrae.fr/hal-04564889

Submitted on 30 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Sizing efficient underdrains for treatment wetlands

Ania Morvannou, Matthieu Dufresne, Marie-Christine Gromaire, Stéphane

Troesch, Nicolas Forquet

► To cite this version:

Ania Morvannou, Matthieu Dufresne, Marie-Christine Gromaire, Stéphane Troesch, Nicolas Forquet. Sizing efficient underdrains for treatment wetlands. Water Science and Technology, 2024, 89 (2), pp.241 - 251. 10.2166/wst.2023.417. hal-04564889

HAL Id: hal-04564889

https://hal.inrae.fr/hal-04564889v1

Submitted on 30 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

IWA
PUBLISHING

# Water Science & Technology

© 2024 The Authors

Water Science & Technology Vol 89 No 2, 241 doi: 10.2166/wst.2023.417

## Sizing efficient underdrains for treatment wetlands

Ania Morvannou IWA *a,*, Matthieu Dufresne<sup>b</sup>, Marie-Christine Gromaire IWA <sup>c</sup>, Stéphane Troesch IWA <sup>d</sup>
and Nicolas Forquet IWA <sup>d</sup>

<figure><img src="image_5.png" /></figure>

AM, 0000-0002-0840-6043; M-CG, 0000-0002-4000-1070; ST, 0000-0002-5353-9187; NF, 0000-0003-1154-5498

## 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 under-sized, 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.

**Key words:** analog model, spatial heterogeneity, stormwater, treatment wetlands, underdrain, urban runoff

## 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.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).

<footnote><sup>a</sup> Rhône, Eco Bird, 3 route du Dôme, Chaponost 69630, France</footnote>

<footnote><sup>b</sup> Bas-Rhin, INSA Strasbourg, 24 Boulevard de la Victoire, Strasbourg 67000, France</footnote>

<footnote><sup>c</sup> Seine-et-Marne, Leesu, Ecole des Ponts, Univ Paris Est Creteil, 6 et 8 av Blaise Pascal, Cité Descartes - Champs sur Marne, Marne-la-Vallee 77455, France</footnote>

<footnote><sup>d</sup> Rhône, INRAE, UR REVERSAAL, 5 rue de la Doua, Villeurbanne 69100, France</footnote>

<footnote>*Corresponding author. E-mail: a.morvannou@ecobird.fr</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

IWA
PUBLISHING

# Water Science & Technology

© 2024 The Authors

Water Science & Technology Vol 89 No 2, 241 doi: 10.2166/wst.2023.417

## Sizing efficient underdrains for treatment wetlands

Ania Morvannou IWA *a,*, Matthieu Dufresne<sup>b</sup>, Marie-Christine Gromaire IWA <sup>c</sup>, Stéphane Troesch IWA <sup>d</sup>
and Nicolas Forquet IWA <sup>d</sup>

<figure><img src="image_4.png" /></figure>

AM, 0000-0002-0840-6043; M-CG, 0000-0002-4000-1070; ST, 0000-0002-5353-9187; NF, 0000-0003-1154-5498

## 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 under-sized, 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.

**Key words:** analog model, spatial heterogeneity, stormwater, treatment wetlands, underdrain, urban runoff

## 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.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).

<footnote><sup>a</sup> Rhône, Eco Bird, 3 route du Dôme, Chaponost 69630, France</footnote>

<footnote><sup>b</sup> Bas-Rhin, INSA Strasbourg, 24 Boulevard de la Victoire, Strasbourg 67000, France</footnote>

<footnote><sup>c</sup> Seine-et-Marne, Leesu, Ecole des Ponts, Univ Paris Est Creteil, 6 et 8 av Blaise Pascal, Cité Descartes - Champs sur Marne, Marne-la-Vallee 77455, France</footnote>

<footnote><sup>d</sup> Rhône, INRAE, UR REVERSAAL, 5 rue de la Doua, Villeurbanne 69100, France</footnote>

<footnote>*Corresponding author. E-mail: a.morvannou@ecobird.fr</footnote>
## GRAPHICAL ABSTRACT

<figure><img src="image_6.png" /></figure>

## 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
## GRAPHICAL ABSTRACT

<figure><img src="image_5.png" /></figure>

## 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 environ-
mental 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 under-
drain 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 perco-
lated flows?

METHODS

Description of the studied site

The system studied in the Life ADSORB project consists in two 600 m² 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.

The underdrain network of TW 1 is made of 160 mm diameter pipes and is distributed as follows (Figure 1): two parallel
antennas with a total length of 105 m merging downstream into a main drainage pipe of 18 m in length. The last 3 m of the
underdrain network are made up of a 315 mm diameter pipe on which the outlet valve is mounted. The slots allowing
the collection of the water along the drainage pipes are oriented upward and are present along the entire length of the
160 mm diameter pipes. They are 10 mm wide, cover one-sixth of the pipe diameter and are spaced every 120 mm
(8 slots/m, 946 slots in total). This represents 84 cm² of opening per meter of pipe.

<figure><img src="image_7.png" /><figcaption>**Figure 1** | 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.</figcaption></figure>
(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 m² 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.

The underdrain network of TW 1 is made of 160 mm diameter pipes and is distributed as follows (Figure 1): two parallel antennas with a total length of 105 m merging downstream into a main drainage pipe of 18 m in length. The last 3 m of the underdrain network are made up of a 315 mm diameter pipe on which the outlet valve is mounted. The slots allowing the collection of the water along the drainage pipes are oriented upward and are present along the entire length of the 160 mm diameter pipes. They are 10 mm wide, cover one-sixth of the pipe diameter and are spaced every 120 mm (8 slots/m, 946 slots in total). This represents 84 cm² of opening per meter of pipe.

<figure><img src="image_6.png" /><figcaption>**Figure 1** | 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.</figcaption></figure>
The annual volume of water to be treated is estimated to be 200,000 m³/year. It is composed of runoff from the Paris ring road (21 ha), permanent parasitic clear water (279 m³/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 m³ 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

First, the linear head losses in the pipes were estimated to evaluate whether the flow through the slots is likely to be evenly distributed. Since the diameter remains unchanged after the junction of the two antennas, significant linear head losses are likely to occur in the main discharge section. It was estimated using the Darcy-Weisbach equation:

$$ \Delta H_L = \lambda \frac{L V^2}{D 2g} \qquad (1) $$

where $\Delta H_L$ is the linear hydraulic head loss in the main drain section (m), $\lambda$ is the Darcy friction coefficient (–), $L$ is the length of the main drain section (m), $D$ is the diameter of the main drain section (m), $V$ is the fluid velocity (m s⁻¹) and $g$ is the gravitational acceleration (m s⁻²).

The Darcy friction coefficient ($\lambda$) is calculated differently depending on the value of the Reynolds number (Re). If Re > 2, 000, it is estimated by Colebrook's equation:

$$ \frac{1}{\sqrt{\lambda}} = 2 \times \log \left[ \frac{k}{3.7D} + \frac{2.51}{Re\sqrt{\lambda}} \right] \qquad (2) $$

where $k$ is the pipe rugosity (m) and $Re$ is the Reynolds number (–). Two values of the rugosity coefficient ($k$) have been chosen: $1.5 \times 10^{-6}$ m corresponding to a new PVC pipe and $1 \times 10^{-4}$ m corresponding to a PVC pipe with some deposit in it (Hager 2010).

If $Re \le 2,000$, only viscosity influences linear hydraulic head losses and $\lambda$ is defined by

$$ \lambda = \frac{64}{Re} \qquad (3) $$

The linear hydraulic head loss in the main drain section is equal to 0.17 m for $k = 1.5 \times 10^{-6}$ m and 0.21 m for $k = 1 \times 10^{-4}$ 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

The singular head loss coefficient of the slot (K) 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 ($P_{upstream}$), from 2,500 to 5,000 Pa for the pressure head downstream of the drain ($P_{downstream}$) and from 0 to $2 \times 10^{-2}$ m³ h⁻¹ for the upstream flow rate ($Q_{in}$). The range of 3,000–7,500 Pa for $P_{upstream}$ was chosen because it represents the amplitude of water pressure variation in the filter. The range of variation for $P_{downstream}$ was chosen based on the $P_{upstream}$ 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
The annual volume of water to be treated is estimated to be 200,000 m³/year. It is composed of runoff from the Paris ring road (21 ha), permanent parasitic clear water (279 m³/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 m³ 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

First, the linear head losses in the pipes were estimated to evaluate whether the flow through the slots is likely to be evenly distributed. Since the diameter remains unchanged after the junction of the two antennas, significant linear head losses are likely to occur in the main discharge section. It was estimated using the Darcy-Weisbach equation:

$$ \Delta H_L = \lambda \frac{L V^2}{D 2g} \qquad (1) $$

where $\Delta H_L$ is the linear hydraulic head loss in the main drain section (m), $\lambda$ is the Darcy friction coefficient (–), $L$ is the length of the main drain section (m), $D$ is the diameter of the main drain section (m), $V$ is the fluid velocity (m s⁻¹) and $g$ is the gravitational acceleration (m s⁻²).

The Darcy friction coefficient ($\lambda$) is calculated differently depending on the value of the Reynolds number (Re). If Re > 2, 000, it is estimated by Colebrook's equation:

$$ \frac{1}{\sqrt{\lambda}} = 2 \times \log \left[ \frac{k}{3.7D} + \frac{2.51}{Re\sqrt{\lambda}} \right] \qquad (2) $$

where $k$ is the pipe rugosity (m) and $Re$ is the Reynolds number (–). Two values of the rugosity coefficient ($k$) have been chosen: 1.5 × 10⁻⁶ m corresponding to a new PVC pipe and 1 × 10⁻⁴ m corresponding to a PVC pipe with some deposit in it (Hager 2010).

If Re ≤ 2,000, only viscosity influences linear hydraulic head losses and $\lambda$ is defined by

$$ \lambda = \frac{64}{Re} \qquad (3) $$

The linear hydraulic head loss in the main drain section is equal to 0.17 m for $k = 1.5 \times 10^{-6}$ m and 0.21 m for $k = 1 \times 10^{-4}$ 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

The singular head loss coefficient of the slot (K) 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 ($P_{upstream}$), from 2,500 to 5,000 Pa for the pressure head downstream of the drain ($P_{downstream}$) and from 0 to $2 \times 10^{-2} m^3 h^{-1}$ for the upstream flow rate ($Q_{in}$). The range of 3,000–7,500 Pa for $P_{upstream}$ was chosen because it represents the amplitude of water pressure variation in the filter. The range of variation for $P_{downstream}$ was chosen based on the $P_{upstream}$ 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
<figure><img src="image_8.png" /><figcaption>**Figure 2** | 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 ($P_{upstream}$) and downstream ($P_{downstream}$) of the drain, and arrow for flow boundary conditions ($Q_{in}$)).</figcaption></figure>

<table><caption>**Table 1** | Boundary conditions and results of CFD simulations</caption><thead><tr><td>P<sub>upstream</sub> (Pa)</td><td>P<sub>downstream</sub> (Pa)</td><td>Pressure head loss (Pa)</td><td>Upstream flow rate Q<sub>in</sub> (m<sup>3</sup> h<sup>-1</sup>)</td><td>Slot flow rate (m<sup>3</sup> h<sup>-1</sup>)</td></tr></thead><tbody><tr><td>3,000</td><td>2,500</td><td>500</td><td>0</td><td>0.19</td></tr><tr><td>3,000</td><td>2,500</td><td>500</td><td>36</td><td>0.18</td></tr><tr><td>3,000</td><td>2,500</td><td>500</td><td>72</td><td>0.16</td></tr><tr><td>7,500</td><td>5,000</td><td>2,500</td><td>0</td><td>0.44</td></tr><tr><td>7,500</td><td>5,000</td><td>2,500</td><td>36</td><td>0.43</td></tr><tr><td>7,500</td><td>5,000</td><td>2,500</td><td>72</td><td>0.41</td></tr><tr><td>7,500</td><td>2,500</td><td>5,000</td><td>0</td><td>0.62</td></tr><tr><td>7,500</td><td>2,500</td><td>5,000</td><td>36</td><td>0.62</td></tr><tr><td>7,500</td><td>2,500</td><td>5,000</td><td>72</td><td>0.60</td></tr></tbody></table>

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 ($K$) is determined by fitting the relationship between the pressure head difference across the slot and the right-hand side of Equation (4).

## 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.

Figure 3 shows the different calculation steps of the analog model. Since flow and pressure are interdependent, a loop is implemented. First, the pressure $P_j$ in the underdrain at slot $j$ is set to a constant value $P_{init}$. 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 $Q(x)$ to the outlet. Knowing the flow rate $Q(x)$, the linear head losses between the slots are then estimated. Finally, $P_j$ are updated based on the linear head losses starting downstream. The loop continues until the differences between the slot level pressures $P_j$ obtained in the previous iteration and those obtained in the current iteration (Delta_P) are less than a given threshold (tolerance).
<figure><img src="image_7.png" /><figcaption>**Figure 2** | 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 ($P_{upstream}$) and downstream ($P_{downstream}$) of the drain, and arrow for flow boundary conditions ($Q_{in}$)).</figcaption></figure>

<table><caption>**Table 1** | Boundary conditions and results of CFD simulations</caption><thead><tr><th>P<sub>upstream</sub> (Pa)</th><th>P<sub>downstream</sub> (Pa)</th><th>Pressure head loss (Pa)</th><th>Upstream flow rate Q<sub>in</sub> (m<sup>3</sup> h<sup>-1</sup>)</th><th>Slot flow rate (m<sup>3</sup> h<sup>-1</sup>)</th></tr></thead><tbody><tr><td>3,000</td><td>2,500</td><td>500</td><td>0</td><td>0.19</td></tr><tr><td>3,000</td><td>2,500</td><td>500</td><td>36</td><td>0.18</td></tr><tr><td>3,000</td><td>2,500</td><td>500</td><td>72</td><td>0.16</td></tr><tr><td>7,500</td><td>5,000</td><td>2,500</td><td>0</td><td>0.44</td></tr><tr><td>7,500</td><td>5,000</td><td>2,500</td><td>36</td><td>0.43</td></tr><tr><td>7,500</td><td>5,000</td><td>2,500</td><td>72</td><td>0.41</td></tr><tr><td>7,500</td><td>2,500</td><td>5,000</td><td>0</td><td>0.62</td></tr><tr><td>7,500</td><td>2,500</td><td>5,000</td><td>36</td><td>0.62</td></tr><tr><td>7,500</td><td>2,500</td><td>5,000</td><td>72</td><td>0.60</td></tr></tbody></table>

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 ($K$) is determined by fitting the relationship between the pressure head difference across the slot and the right-hand side of Equation (4).

## 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.

Figure 3 shows the different calculation steps of the analog model. Since flow and pressure are interdependent, a loop is implemented. First, the pressure $P_j$ in the underdrain at slot $j$ is set to a constant value $P_{init}$. 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 $Q(x)$ to the outlet. Knowing the flow rate $Q(x)$, the linear head losses between the slots are then estimated. Finally, $P_j$ are updated based on the linear head losses starting downstream. The loop continues until the differences between the slot level pressures $P_j$ obtained in the previous iteration and those obtained in the current iteration (Delta_P) are less than a given threshold (tolerance).
<figure><img src="image_9.png" /><figcaption>**Figure 3** | Analog model flow chart.</figcaption></figure>

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):

$$ \Delta H_s = K \frac{Q_s^2}{2gS_s^2} \qquad (4) $$

where $\Delta H_s$ is the singular hydraulic head loss (m), $K$ is the head loss coefficient of the slot (–), $Q_s$ is the flow rate passing through the slot ($m^3 s^{-1}$) and $S_s$ is the cross-sectional area of the slot ($m^2$). The estimation of the head loss coefficient ($K$) 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 ($H_{outlet}$) is calculated by estimating the singular head losses associated to the cumulative flow rate passing through the throttle outflow,

$$ Q_{outlet} = C_d S \sqrt{2g(H_{outlet} - H_{throttle})} \qquad (5) $$
<figure><img src="image_8.png" /><figcaption>**Figure 3** | Analog model flow chart.</figcaption></figure>

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):

$$ \Delta H_s = K \frac{Q_s^2}{2gS_s^2} \qquad (4) $$

where $\Delta H_s$ is the singular hydraulic head loss (m), $K$ is the head loss coefficient of the slot (–), $Q_s$ is the flow rate passing through the slot ($m^3 s^{-1}$) and $S_s$ is the cross-sectional area of the slot ($m^2$). The estimation of the head loss coefficient ($K$) 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 ($H_{outlet}$) is calculated by estimating the singular head losses associated to the cumulative flow rate passing through the throttle outflow,

$$ Q_{outlet} = C_d S \sqrt{2g(H_{outlet} - H_{throttle})} \qquad (5) $$
where $Q_{outlet}$ is the flow rate exiting the drain (calculated by summing up the incoming flow rate passing through each slot) (m³ s⁻¹), $C_d$ is the singular head loss coefficient (equal to 0.62 for a circular opening) (–) and $H_{throttle}$ is the height of the throttle opening measured from the bottom of the TW (m).

* 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 equation

$$H_i = (1 - \alpha)H_{i-1} + \alpha H_{i-2} \quad (6)$$

where $\alpha$ 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 i – 1 and the loop can be exit when

$$\max |H_i - H_{i-1}| \le tol \quad (7)$$

where tol is the tolerance, set to $1 \times 10^{-3}$ in this study.

## Sizing optimization algorithm

Finally, with the appropriate singular head loss coefficient (K) 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 ($f_{obj}$) has been defined that measures the difference in flow rate between the first and the last slot of each branch,

$$f_{obj} = \frac{|Q_{first} - Q_{last}|}{Q_{last}}$$

with $Q_{first}$ being the flow rate passing through the first slot of a branch (m³ s⁻¹) and $Q_{last}$ being the flow rate passing through the last slot of a branch (m³ s⁻¹).

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 ($2.5 \times 10^{-1}$ 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

Figure 4 displays the results of the slot numerical model, calculating the head loss across the slot versus the flow crossing it to the power of two for all combinations of boundary conditions tested (see Table 1). They exhibit a linear trend whose slope corresponds to the singular head loss of the slot: $K = 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).

### Estimation of flow distribution along the underdrain network using the analog model

To illustrate the influence of linear hydraulic head loss on the flow distribution across the slots of the underdrain network, the analog model was run for water heights above the underdrain network ($H_{filter}$) 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
where $Q_{outlet}$ is the flow rate exiting the drain (calculated by summing up the incoming flow rate passing through each slot) (m³ s⁻¹), $C_d$ is the singular head loss coefficient (equal to 0.62 for a circular opening) (–) and $H_{throttle}$ is the height of the throttle opening measured from the bottom of the TW (m).

* 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 equation

$$H_i = (1 - \alpha)H_{i-1} + \alpha H_{i-2} \quad (6)$$

where $\alpha$ 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 i – 1 and the loop can be exit when

$$\max |H_i - H_{i-1}| \le tol \quad (7)$$

where tol is the tolerance, set to $1 \times 10^{-3}$ in this study.

## Sizing optimization algorithm

Finally, with the appropriate singular head loss coefficient (K) 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 ($f_{obj}$) has been defined that measures the difference in flow rate between the first and the last slot of each branch,

$$f_{obj} = \frac{|Q_{first} - Q_{last}|}{Q_{last}}$$

with $Q_{first}$ being the flow rate passing through the first slot of a branch (m³ s⁻¹) and $Q_{last}$ being the flow rate passing through the last slot of a branch (m³ s⁻¹).

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 ($2.5 \times 10^{-1}$ 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

Figure 4 displays the results of the slot numerical model, calculating the head loss across the slot versus the flow crossing it to the power of two for all combinations of boundary conditions tested (see Table 1). They exhibit a linear trend whose slope corresponds to the singular head loss of the slot: $K = 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).

### Estimation of flow distribution along the underdrain network using the analog model

To illustrate the influence of linear hydraulic head loss on the flow distribution across the slots of the underdrain network, the analog model was run for water heights above the underdrain network ($H_{filter}$) 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
<figure><img src="image_10.png" /><figcaption>**Figure 4** | Determination of the singular head loss of the slot by fitting results of numerical experiments.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 5** | 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 $k = 1.5 \times 10^{-6}$ m).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 4** | Determination of the singular head loss of the slot by fitting results of numerical experiments.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5** | 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 $k = 1.5 \times 10^{-6}$ m).</figcaption></figure>
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.

## Optimized sizing of the underdrain network

In order to reduce the linear head losses, pipe diameters were optimized using the algorithm described in “Sizing optimization algorithm” paragraph in the “Methods” section. Calculations were performed for pipes with smooth and rough rugosity ($k = 1.5 \times 10^{-6} m$ and $k = 1 \times 10^{-4} m$, 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.

Finally, Figure 6 is the result of applying the analog model to the optimized underdrain network with a rugosity value of $k = 1.5 \times 10^{-6} m$. 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.

<table><caption>**Table 2** | Drainage pipe diameters in mm: current values versus optimized ones for smooth and rough pipes</caption><thead><tr><td>Pipe</td><td>Current diameter (mm)</td><td>Optimized diameter (mm) (k = 1.5 × 10<sup>-6</sup> m)</td><td>Optimized diameter (mm) (k = 1 × 10<sup>-4</sup> m)</td></tr></thead><tbody><tr><td>Main</td><td>160</td><td>315</td><td>315</td></tr><tr><td>Antennas</td><td>160</td><td>250</td><td>250</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 6** | 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 $k = 1.5 \times 10^{-6} m$).</figcaption></figure>
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.

## Optimized sizing of the underdrain network

In order to reduce the linear head losses, pipe diameters were optimized using the algorithm described in “Sizing optimization algorithm” paragraph in the “Methods” section. Calculations were performed for pipes with smooth and rough rugosity ($k = 1.5 \times 10^{-6} m$ and $k = 1 \times 10^{-4} m$, 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.

Finally, Figure 6 is the result of applying the analog model to the optimized underdrain network with a rugosity value of $k = 1.5 \times 10^{-6} m$. 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.

<table><caption>**Table 2** | Drainage pipe diameters in mm: current values versus optimized ones for smooth and rough pipes</caption><thead><tr><td>Pipe</td><td>Current diameter (mm)</td><td>Optimized diameter (mm) ($k = 1.5 \times 10^{-6} m$)</td><td>Optimized diameter (mm) ($k = 1 \times 10^{-4} m$)</td></tr></thead><tbody><tr><td>Main</td><td>160</td><td>315</td><td>315</td></tr><tr><td>Antennas</td><td>160</td><td>250</td><td>250</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 6** | 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 $k = 1.5 \times 10^{-6} m$).</figcaption></figure>
## 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.

## REFERENCES

Adyel, T. M., Oldham, C. E. & Hipsey, M. R. 2016 *Stormwater nutrient attenuation in a constructed wetland with alternating surface and subsurface flow pathways: Event to annual dynamics*. *Water Research* **107**, 66–82. https://doi.org/10.1016/j.watres.2016.10.005.

Biswal, B. K., Bolan, N., Zhu, Y.-G. & Balasubramanian, R. 2022 *Nature-based systems (NbS) for mitigation of stormwater and air pollution in urban areas: A review*. *Resources, Conservation and Recycling* **186**, 106578. https://doi.org/10.1016/j.resconrec.2022.106578.

Choi, J. Y., Maniquiz, M. C., Geronimo, F. K., Lee, S. Y., Lee, B. S. & Kim, L. H. 2012 *Development of a horizontal subsurface flow modular constructed wetland for urban runoff treatment*. *Water Science and Technology* **66** (9), 1950–1957. https://doi.org/10.2166/wst.2012.306.

DWA-A 178. 2017 *Arbeitsblatt DWA-A 178: Empfehlungen Für Planung, Bau Und Betrieb von Retentionsbodenfiltern Zur Weitergehenden Regenwasserbehandlung Im Misch- Und Trennsystem (Recommendations for Planning, Construction and Operation of Retention Soil Filters for Advanced Stormwater Treatment in Combined and Separate Sewer Systems)*. German Association for Water Management, Wastewater; Waste, GFA, Hennef (in German; draft).

Forquet, N. & Dufresne, M. 2015 *Simple deterministic model of the hydraulic buffer effect in septic tanks*. *Water and Environment Journal* **29** (3), 360–364. https://doi.org/10.1111/wej.12114.

Hager, W. H. 2010 *Wastewater Hydraulics*. Springer, Berlin. https://doi.org/10.1007/978-3-642-11383-3.
## 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.

## REFERENCES

Adyel, T. M., Oldham, C. E. & Hipsey, M. R. 2016 *Stormwater nutrient attenuation in a constructed wetland with alternating surface and subsurface flow pathways: Event to annual dynamics*. *Water Research* **107**, 66–82. https://doi.org/10.1016/j.watres.2016.10.005.

Biswal, B. K., Bolan, N., Zhu, Y.-G. & Balasubramanian, R. 2022 *Nature-based systems (NbS) for mitigation of stormwater and air pollution in urban areas: A review*. *Resources, Conservation and Recycling* **186**, 106578. https://doi.org/10.1016/j.resconrec.2022.106578.

Choi, J. Y., Maniquiz, M. C., Geronimo, F. K., Lee, S. Y., Lee, B. S. & Kim, L. H. 2012 *Development of a horizontal subsurface flow modular constructed wetland for urban runoff treatment*. *Water Science and Technology* **66** (9), 1950–1957. https://doi.org/10.2166/wst.2012.306.

DWA-A 178. 2017 *Arbeitsblatt DWA-A 178: Empfehlungen Für Planung, Bau Und Betrieb von Retentionsbodenfiltern Zur Weitergehenden Regenwasserbehandlung Im Misch- Und Trennsystem (Recommendations for Planning, Construction and Operation of Retention Soil Filters for Advanced Stormwater Treatment in Combined and Separate Sewer Systems)*. German Association for Water Management, Wastewater; Waste, GFA, Hennef (in German; draft).

Forquet, N. & Dufresne, M. 2015 *Simple deterministic model of the hydraulic buffer effect in septic tanks*. *Water and Environment Journal* **29** (3), 360–364. https://doi.org/10.1111/wej.12114.

Hager, W. H. 2010 *Wastewater Hydraulics*. Springer, Berlin. https://doi.org/10.1007/978-3-642-11383-3.
Liu, J., Sample, D. J., Bell, C. & Guan, Y. 2014 Review and research needs of bioretention used for the treatment of urban stormwater. *Water* **6** (4), 1069–1099. https://doi.org/10.3390/w6041069.

Meyer, D., Molle, P., Esser, D., Troesch, S., Masi, F. & Dittmer, U. 2013 Constructed wetlands for combined sewer overflow treatment – Comparison of German, French and Italian approaches. *Water* **5** (1), 1–12. https://doi.org/10.3390/w5010001.

Molle, P., Fournel, J., Meyer, D., Troesch, S., Clement, F., Brelot, E., Bacot, L., Guillermond, S., de Brito, C., Toussaint, J.Y., Vareilles, S., Ah Leung, S., Lipeme Kouyi, G., Bichet, Q., Chocat, B. & Esser, D. 2013 *Systèmes Extensifs Pour La Gestion Et Le Traitement Des Eaux Urbaines de Temps de Pluie* (Extensive Systems for Storm Water and Combine Sewer Overflow Treatment). Agence Nationale de la Recherche. Available from: https://hal.inrae.fr/hal-02599141v1.

Pálfy, T. G., Gerodolle, M., Gourdon, R., Meyer, D., Troesch, S. & Molle, P. 2017 Performance assessment of a vertical flow constructed wetland treating unsettled combined sewer overflow. *Water Science and Technology* **75** (11), 2586–2597. https://doi.org/10.2166/wst.2017.126.

Rizzo, A., Tondera, K., Pálfy, T. G., Dittmer, U., Meyer, D., Schreiber, C., Zacharias, N., Ruppelt, J.P., Esser, D., Molle, P., Troesch, S. & Masi, F. 2020 Constructed wetlands for combined sewer overflow treatment: A State-of-the-Art review. *Science of the Total Environment* **727**, 138618. https://doi.org/10.1016/j.scitotenv.2020.138618.

Roy-Poirier, A., Champagne, P. & Filion, Y. 2010 Review of bioretention system research and design: Past, present, and future. *Journal of Environmental Engineering* **136** (9), 878–889. https://doi.org/10.1061/(ASCE)EE.1943-7870.0000227.

Schmitt, N., Wanko, A., Laurent, J., Bois, P., Molle, P. & Mosé, R. 2015 Constructed wetlands treating stormwater from separate sewer networks in a residential Strasbourg urban catchment area: Micropollutant removal and fate. *Journal of Environmental Chemical Engineering* **3** (4, Part A), 2816–2824. https://doi.org/10.1016/j.jece.2015.10.008.

Shutes, R. B. E., Revitt, D. M., Lagerberg, I. M. & Barraud, V. C. E. 1999 The design of vegetative constructed wetlands for the treatment of highway runoff. *Science of the Total Environment* **235** (1), 189–197. https://doi.org/10.1016/S0048-9697(99)00212-0.

Tondera, K. 2019 Evaluating the performance of constructed wetlands for the treatment of combined sewer overflows. *Ecological Engineering* **137**, 53–59. https://doi.org/10.1016/j.ecoleng.2017.10.009.

Vymazal, J., Kröpfelová, L., 2008 *Wastewater Treatment in Constructed Wetlands with Horizontal Sub-Surface Flow*, Vol. 14. (Alloway, B. J. & Trevors, J. T., eds). Environmental Pollution. Springer Netherlands, Dordrecht. https://doi.org/10.1007/978-1-4020-8580-2.

Walaszek, M., Bois, P., Laurent, J., Lenormand, E. & Wanko, A. 2018 Urban stormwater treatment by a constructed wetland: Seasonality impacts on hydraulic efficiency, physico-chemical behavior and heavy metal occurrence. *Science of the Total Environment* **637–638**, 443–454. https://doi.org/10.1016/j.scitotenv.2018.04.325.

First received 30 June 2023; accepted in revised form 1 December 2023. Available online 22 December 2023
Liu, J., Sample, D. J., Bell, C. & Guan, Y. 2014 Review and research needs of bioretention used for the treatment of urban stormwater. *Water* **6** (4), 1069–1099. https://doi.org/10.3390/w6041069.

Meyer, D., Molle, P., Esser, D., Troesch, S., Masi, F. & Dittmer, U. 2013 Constructed wetlands for combined sewer overflow treatment – Comparison of German, French and Italian approaches. *Water* **5** (1), 1–12. https://doi.org/10.3390/w5010001.

Molle, P., Fournel, J., Meyer, D., Troesch, S., Clement, F., Brelot, E., Bacot, L., Guillermard, S., de Brito, C., Toussaint, J.Y., Vareilles, S., Ah Leung, S., Lipeme Kouyi, G., Bichet, Q., Chocat, B. & Esser, D. 2013 *Systèmes Extensifs Pour La Gestion Et Le Traitement Des Eaux Urbaines de Temps de Pluie* (Extensive Systems for Storm Water and Combine Sewer Overflow Treatment). Agence Nationale de la Recherche. Available from: https://hal.inrae.fr/hal-02599141v1.

Pálfy, T. G., Gerodolle, M., Gourdon, R., Meyer, D., Troesch, S. & Molle, P. 2017 Performance assessment of a vertical flow constructed wetland treating unsettled combined sewer overflow. *Water Science and Technology* **75** (11), 2586–2597. https://doi.org/10.2166/wst.2017.126.

Rizzo, A., Tondera, K., Pálfy, T. G., Dittmer, U., Meyer, D., Schreiber, C., Zacharias, N., Ruppelt, J.P., Esser, D., Molle, P., Troesch, S. & Masi, F. 2020 Constructed wetlands for combined sewer overflow treatment: A State-of-the-Art review. *Science of the Total Environment* **727**, 138618. https://doi.org/10.1016/j.scitotenv.2020.138618.

Roy-Poirier, A., Champagne, P. & Filion, Y. 2010 Review of bioretention system research and design: Past, present, and future. *Journal of Environmental Engineering* **136** (9), 878–889. https://doi.org/10.1061/(ASCE)EE.1943-7870.0000227.

Schmitt, N., Wanko, A., Laurent, J., Bois, P., Molle, P. & Mosé, R. 2015 Constructed wetlands treating stormwater from separate sewer networks in a residential Strasbourg urban catchment area: Micropollutant removal and fate. *Journal of Environmental Chemical Engineering* **3** (4, Part A), 2816–2824. https://doi.org/10.1016/j.jece.2015.10.008.

Shutes, R. B. E., Revitt, D. M., Lagerberg, I. M. & Barraud, V. C. E. 1999 The design of vegetative constructed wetlands for the treatment of highway runoff. *Science of the Total Environment* **235** (1), 189–197. https://doi.org/10.1016/S0048-9697(99)00212-0.

Tondera, K. 2019 Evaluating the performance of constructed wetlands for the treatment of combined sewer overflows. *Ecological Engineering* **137**, 53–59. https://doi.org/10.1016/j.ecoleng.2017.10.009.

Vymazal, J., Kröpfelová, L., 2008 *Wastewater Treatment in Constructed Wetlands with Horizontal Sub-Surface Flow*, Vol. 14. (Alloway, B. J. & Trevors, J. T., eds). Environmental Pollution. Springer Netherlands, Dordrecht. https://doi.org/10.1007/978-1-4020-8580-2.

Walaszek, M., Bois, P., Laurent, J., Lenormand, E. & Wanko, A. 2018 Urban stormwater treatment by a constructed wetland: Seasonality impacts on hydraulic efficiency, physico-chemical behavior and heavy metal occurrence. *Science of the Total Environment* **637–638**, 443–454. https://doi.org/10.1016/j.scitotenv.2018.04.325.

First received 30 June 2023; accepted in revised form 1 December 2023. Available online 22 December 2023