# Regression-based design optimization of French treatment wetlands

Zoé Legeai, Pascal Molle, Alberto Tonda, Nicolas Forquet

► To cite this version:

Zoé Legeai, Pascal Molle, Alberto Tonda, Nicolas Forquet. Regression-based design optimization of French treatment wetlands. Water Science and Technology, 2025, 91 (11), pp.1248 - 1264. 10.2166/wst.2025.071. hal-05226448

HAL Id: hal-05226448

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

Submitted on 27 Aug 2025

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_2.png" /></figure>

Water Science & Technology

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

IWA
PUBLISHING

© 2025 The Authors

Water Science & Technology Vol 91 No 11, 1248 doi: 10.2166/wst.2025.071

Regression-based design optimization of French treatment wetlands

Zoé Legeai<sup>a,†</sup>, Pascal Molle<sup>IWA</sup> <sup>a,†</sup>, Alberto Tonda<sup>b,c</sup> and Nicolas Forquet<sup>IWA</sup> <sup>a,∗</sup>

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

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

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

ZL, 0009-0003-2550-8206; PM, 0000-0003-3364-6118; AT, 0000-0001-5895-4809; NF, 0000-0003-1154-5498

ABSTRACT

This study presents a novel approach for the automatic optimization of sizing French vertical flow wetlands, a widely used solution in waste-water treatment across France, with over 6,000 operational units. By leveraging a comprehensive database collected by the French National Research Institute for Agriculture, Food and Environment (INRAE), we have developed robust regression models that link total Kjeldahl nitrogen to surface and total chemical oxygen demand to material depth. By integrating these models into a unified optimization framework, our methodology enhances the design of the entire treatment train, particularly addressing how the sizing of the first stage influences the second stage sizing while complying with stringent effluent quality standards. Extensive validation confirms the effectiveness and reliability of our approach, which is tailored to specific influent characteristics and regulatory requirements. Our strategy improves upon traditional sizing methods, which often rely on surface loads, by allowing for refined optimization based on specific performance targets. The findings contribute to improved design practices and suggest potential directions for future developments, including the integration of additional treatment processes and multi-criteria decision analysis for broader applicability in wastewater management.

**Key words:** COD and TKN removal, regression models, treatment train optimization, uncertainty analysis, vertical flow treatment wetlands, wastewater treatment

HIGHLIGHTS

* A novel optimization framework for sizing treatment wetlands considering the entire treatment train.

* It combines regression models with global optimization under constraints.

* Uncertainties are provided for a more secure design.

* Filter depths over 30 cm are necessary for stringent outlet levels or heavily loaded influents.

* Sizing for TKN outlet levels down to 10 mg/L is achievable with acceptable uncertainty.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

<footnote><sup>a</sup> INRAE, REVERSAAL Research Unit, Villeurbanne, France</footnote>

<footnote><sup>b</sup> INRAE, MIA-PS Research Unit, Université Paris-Saclay, Palaiseau, France</footnote>

<footnote><sup>c</sup> CNRS, Institut des Systèmes Complexes Paris Ile-de-France, Paris, France</footnote>

<footnote>*Corresponding author. E-mail: nicolas.forquet@inrae.fr</footnote>
## GRAPHICAL ABSTRACT

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

## INTRODUCTION

The French vertical flow wetland (Dotro et al. 2021) is a two-stage design whose main feature is the direct discharge of raw wastewater in its first treatment stage (Molle et al. 2005). The first experiments in France were initiated in the 1980s, and the design rules were derived empirically by observing the successive developments (Liénard 2010; Morvannou et al. 2015).

At present, there are over 6,000 operational French vertical flow wetlands serving small communities in France and a growing number abroad (Ministère de la Transition Écologique 2023). The regulatory monitoring results have been collected in a database (Morvannou et al. 2015) by the French National Research Institute for Agriculture, Food and Environment (INRAE). This includes the collection of accurate design data for approximately 100 wastewater treatment plants, which has been used to define design equations based on regression models.

Further data collection has enhanced existing data-driven models, particularly through the incorporation of filter depth (Millot et al. 2016). Furthermore, the advent of more stringent outflow level quality standards has also prompted the necessity for these improvements. Regression models represent one among several other numerical models that have been developed over the last three decades with the objective of understanding and predicting the performance of treatment wetlands (Meyer et al. 2015; Defo et al. 2017; Yuan et al. 2020).

These models can be classified into two categories:

- Data-driven models. These models are typically concerned with forecasting outcomes rather than processes, focusing on inputs and outputs based solely on the collected data. These models range from linear regressions to simple machine learning regressors to deep learning methods (Nayak et al. 2006; Naz et al. 2009).

- Mechanistic process-based models. These models are designed to provide a more detailed description of the various mechanisms of contaminant transformation and degradation that occur at the core of a treatment wetland. The hydrodynamic is represented either using tank-in-series models (Kadlec & Knight 1996; Kadlec 2000; Stone et al. 2004; Stein et al. 2007; Sun et al. 2008; Kadlec & Wallace 2009) or by solving the Richards equation for water flow in variably saturated
porous media, using for example the finite element method (Langergraber & Šimůnek 2005). Biodegradation is estimated using either first-order models (Rousseau et al. 2004) or Monod models (Mitchell & McNevin 2001). Various models and software packages are available for use, with the HYDRUS wetland module (Langergraber & Šimůnek 2005, 2012; García et al. 2010; Morvannou et al. 2014) being the most widely used. Other favored choices include FITOVERT (Brovelli et al. 2007; Giraldi et al. 2010), BIO_PORE (Samsó & Garcia 2013a, b), among others.

The majority of model applications have concentrated on predicting performance, with only a small number of studies focusing on design optimization (Pálfy et al. 2016; Ruiz 2017). Furthermore, the studies that have addressed this issue have focused on a specific type of wetland. However, the majority of wastewater treatment plants consist of multiple treatment wetlands arranged in a train. When considering the optimization of trains, it is imperative to consider the issue of combining reactors, which is often overlooked and can present challenges in terms of dimensionality. Optimization strategies for comprehensive treatment trains have been simultaneously investigated for conventional treatment systems (Grau et al. 2007; Rieger et al. 2010).

The initial hurdle in optimizing wetland design is the lack of high-quality performance data (Guillaume-Ruty et al. 2024), restricting the use of both empirical and mechanistic models. This research tackles this issue by enhancing current empirical models with new data derived from French vertical flow wetlands. Furthermore, unlike earlier studies that usually examine each treatment stage separately, our methodology considers the entire treatment train, including the interactions between the first and second stages. Through the integration of robust regression models into an optimization algorithm, we present an automated strategy to identify the minimum reactor volume necessary to meet the desired effluent quality standards.

This study has two main objectives:

1. To develop an optimization framework that refines the sizing of French vertical flow wetlands using regression-based models.

2. To assess the impact of uncertainty on design decisions, improving the reliability of optimized solutions.

This study thereby enhances the design practices for wetlands as well as the general domain of wastewater treatment optimization, offering potential applications that extend beyond French vertical flow wetlands.

# MATERIALS AND METHODS

## Regression model

The proposed design support model utilizes regression models implemented using the non-linear least squares (nls) function in R (R Core Team 2024) for each treatment stage. Each pollutant's removal was modeled independently for each stage, with the complexity of the models varying according to the pollutants of interest, namely total suspended solids (TSS), biochemical oxygen demand (BOD5), total chemical oxygen demand (CODt), and total Kjeldahl nitrogen (TKN).

The regression correlations were established following the works of Molle et al. (2005, 2008), with some adaptations for CODt and TKN.

## Total suspended solids

Regarding TSS, we applied the equation from Molle et al. (2005), which assumes a linear relationship between the inlet and outlet load. The reduction equations are as follows:

$$M_{\text{out 1}} = M_{\text{in 1}} \times 0.1 \qquad (1a)$$

$$M_{\text{out 2}} = M_{\text{in 2}} \times 0.2 \qquad (1b)$$

where $M_{in i}$ is the mass load at the inlet of the stage, here for TSS (gTSS/m²/d); $M_{out i}$ is the mass load at the outlet of the stage, here for TSS (gTSS/m²/d); and i is the stage number.
## Biochemical oxygen demand

For BOD₅, similar to TSS, we employed the equation from Molle et al. (2005). The reduction equations are as follows:

$$M_{\text{out 1}} = M_{\text{in 1}} \times 0.1 \qquad (2a)$$

$$M_{\text{out 2}} = M_{\text{in 2}} \times 0.25 \qquad (2b)$$

where $M_{\text{in i}}$ is the mass load at the inlet of the stage, here for BOD₅ (gBOD₅/m²/d); $M_{\text{out i}}$ is the mass load at the outlet of the stage, here for BOD₅ (gBOD₅/m²/d); and i is the stage number.

## Total Kjeldahl nitrogen

While the removals of CODₜ, BOD₅, and TSS demonstrate a linear relationship with inlet mass loads ($M_{\text{in}}$), this is not the case for TKN whose treatment efficiency diminishes as influent loads increase (Molle et al. 2005, 2008). As a result, the TKN treatment objective has a significant impact on the determination of filter surface area.

The data from Molle et al. (2005, 2008) were refitted using a non-linear fitting technique which was preferred over log-transformation (Tedoldi et al. 2025) in order to establish the removal equation. The resulting equations are as follows:

$$M_{\text{out 1}} = 1.1375 \times M_{\text{in 1}}^{0.8168} \qquad (3a)$$

$$M_{\text{out 2}} = 1.5093 \times M_{\text{in 2}}^{0.7887} \qquad (3b)$$

where $M_{\text{in i}}$ is the mass load at the inlet of the stage, here for TKN (gTKN/m²/d); $M_{\text{out i}}$ is the mass load at the outlet of the stage, here for TKN (gTKN/m²/d); and i is the stage number.

These are analogous to the equations presented by Molle et al. (2005, 2008), but offer a superior fit and are less prone to bias in the upper range of TKN loads or when very stringent treatment objectives are required. Please refer to Table S1 in Supplementary Material 1 for a detailed display of the original literature equations, corrected versions, and R² coefficient of determination. For further insight, an uncertainty analysis to visualize the propagation of regression errors is also available in Supplementary Material 1, Figure S1.

## Total chemical oxygen demand

For CODₜ, the standard design approach assumes that removal performance is stable and independent of filter depth. However, recent research has demonstrated that filter depth can significantly impact dissolved biodegradable COD (COD₉₉) removal in first stage filters (Millot et al. 2016) and total organic carbon (TOC) removal in sand-based vertical filters (Olsson 2012). The proposed model incorporates the impact of depth on COD₉₉ removal, thereby enabling the calculation of the total volume of the filter. This model integrates findings from Millot et al. (2016) and Olsson (2012) and assumes that COD₉₉ undergoes first-order degradation as a function of depth.

According to Olsson (2012), dissolved TOC values were converted into dissolved COD (CODₐ). Several studies have previously established relationships between TOC and CODₜ concentrations (Meijer 2004; Henze et al. 2008; Dubber & Gray 2010). The equation developed by Meijer (2004) was selected for its simplicity and its alignment with the upper limit of the range of the CODₜ/TOC ratio as reported by Henze et al. (2008) ([CODₜ]/[TOC] ranges from 2 to 2.5 gCODₜ/gTOC). The resulting ratio is expressed as follows: [TOC]/[CODₜ] = 0.4 gTOC/gCODₜ.

As previously mentioned, these measurements refer to the COD₉₉. To ensure consistency with the actual measurements of CODₜ, a CODₜ fractionation was performed between the COD₉₉, dissolved inert COD (CODₐ), and particulate COD (CODₚ). The CODₐ was assumed to represent 4% of the incoming CODₜ concentration (Gillot & Choubert 2010), with a maximum limit set at 30 mgCODₐ/L. The particulate fraction was calculated based on the inlet concentration of TSS according to the following formula (INRAE database):

$$[\text{COD}_p]_{\text{in 1}} = [\text{TSS}]_{\text{in 1}} \times 1.1 \qquad (4)$$

where $[\text{COD}_p]_{\text{in 1}}$ is the concentration of the particulate COD at the entrance of the first stage (mgCODₚ/L); and $[\text{TSS}]_{\text{in 1}}$ is the concentration of TSS at the entrance of the first stage (mgTSS/L).

The equation used to model the removal of the CODₚ is identical to that used for TSS.
The particulate and inert fractions of CODt were incorporated into the regression by combining them into a single [COD]* term. This step allows for the estimation of the removal rate of only the dissolved biodegradable fraction. Accordingly, 30 mgCOD*/L was deducted from the work of Millot et al. (2016) for the inert biodegradable fraction, while 23 mgCOD*/L was subtracted from Olsson (2012). This adjustment ensured consistency in the ratio of [COD]*/[CODdb] compared with the work of Millot et al. (2016). The experimental data were then fitted using the following model:

$$[COD_{db}]_{out} = [COD_{db}]_{in} \times e^{a \times Z_{i}} \qquad (5)$$

where [CODdb]out i is the concentration of the dissolved biodegradable COD fraction at the outlet of the stage (mgCODdb/L); [CODdb]in i is the concentration of the dissolved biodegradable COD fraction at the inlet of the stage (mgCODdb/L); a is the coefficient to be determined by non-linear regression; Zi is the depth of the first material layer (m); and i is the stage number.

To verify the accuracy of the regression model, we used the R² statistical metric. Additionally, leave-one-out (LOO) cross-validation was performed to estimate the standard deviations for each stage (for further details, see James et al. 2013).

All equations used for the different pollutants are summarized in Figure 1 (simplified version; the complete version is available in Supplementary Material 2, Figure S2).

<figure><img src="image_8.png" /><figcaption>**Figure 1** | Summary of the various equations and numerical values used in our optimization methodology (simplified version). The abbreviations used for the parameters correspond to the same elements presented in the equations above. The complete version is available in Supplementary Material 2, Figure S2.</figcaption></figure>
## Validation

The model was validated using the INRAE database (Morvannou et al. 2015) for each stage individually. The database contains data from approximately 700 treatment wetlands with various configurations. Of these, 67% are two-stage French vertical flow wetlands arranged in series, while 6% consist of only one stage of French vertical flow wetlands. A sample of 415 French vertical flow wetlands plants was selected to assess treatment performance, using 24-h flow-proportional composite samples collected from both the inlet and outlet. These samples included measurements of COD, BOD, TSS, and TKN, conducted according to APHA (2012) standards. The data from the 24-h composite samples were verified against standard literature ranges and typical ratios utilizing the Chauvenet criterion (Guillaume-Ruty et al. 2024). Samples lacking flow rate measurements at the inlet or outlet were excluded.

## Constraints and bounds

To propose a coherent and relevant optimization problem, it is necessary to first establish sizing rules expressed as constraints or bounds. There are three distinct categories: daily surface load constraints, outlet concentration constraints, and technical bounds.

The distinction between a constraint and a bound depends entirely on the optimization algorithm. Bounds prevent the optimization algorithm from generating candidate solutions with values outside of the specified limits. Constraints are conditions that impose a severe penalty on the quality of the candidate solution, but in principle solutions, violating constraints can be produced during the optimization process, even if the correct solution should ideally abide by all given constraints.

### Daily surface load constraints

Dotro et al. (2021) provided recommendations on the maximum pollutant load allowable in French vertical flow wetlands. These recommendations are listed in Table 1.

Note that the values for TKN differ slightly from those recommended by Dotro et al. (2021), who suggest 30 gTKN/m²/d for the first stage and 15 gTKN/m²/d for the second stage. The TKN load constraints indicated here are less stringent than those indicated by Dotro et al. (2021). The previous approach was to ensure the reliability of the system by setting very stringent values. However, as we are optimizing, we can set less stringent values for cases where low nitrification is required. For the remaining parameters (TSS, BOD, and COD), we align with established literature values to mitigate the risk of long-term clogging.

The inlet load constraints are expressed as:

$$ \text{Maximum load}_i \ge Q \times \frac{[\text{pollutant}]_{\text{in } i}}{S_i} \qquad (6) $$

where Maximum load, $i$ is the maximal load of the pollutant of concern for the stage (g/m²/d); $Q$ is the flow rate (m³/d); $[\text{pollutant}]_{\text{in } i}$ is the inlet concentration of the pollutant of concern for the stage (mg/L); $S_i$ is the surface area of the stage (m²); and $i$ is the stage number.

### Outlet concentration constraint

Regarding treated effluent limits, the legal regulations in many countries focus solely on the concentration of pollutants. Consequently, constraints are established based on this parameter. The outlet concentrations are calculated using the various removal equations presented above. These values are contingent on the surface area and depth of each stage, as well as inlet concentrations and flow rate. It is imperative that these values fall below objective concentration values stipulated by the user. This inequality must be verified for each of the pollutants studied (TSS, BOD, TKN, and COD). The following

<table><caption>**Table 1** | Daily surface loads for classical French vertical flow wetlands design under dry weather conditions</caption><thead><tr><td>Treatment stage</td><td>TSS (gTSS/m²/d)</td><td>BOD<sub>5</sub> (gBOD<sub>5</sub>/m²/d)</td><td>TKN (gTKN/m²/d)</td><td>COD<sub>t</sub> (gCOD<sub>t</sub>/m²/d)</td></tr></thead><tbody><tr><td>First stage</td><td>150</td><td>150</td><td>50</td><td>350</td></tr><tr><td>Second stage</td><td>30</td><td>20</td><td>25</td><td>70</td></tr></tbody></table>

Values given are per square meter of bed in operation.
equation can therefore be established:

$$ [\text{pollutant}]_{\text{out } 2} = \begin{cases} f(S_1, Z_1, S_2, Z_2, [\text{pollutant}]_{\text{in } 1}, Q) \\ \le [\text{pollutant}]_{\text{obj}} \end{cases} \quad (7) $$

where $S_i$ is the surface area of the stage (m²); $Z_i$ is the depth of the stage (m²); $[\text{pollutant}]_{\text{out } 2}$ is the outlet concentration of the pollutant of concern for the treatment train (mg/L); $[\text{pollutant}]_{\text{in } 1}$ is the inlet concentration of the pollutant of concern for the treatment train (mg/L); $[\text{pollutant}]_{\text{obj}}$ is the objective concentration of the pollutant of concern for the treatment train (mg/L); Q is the flow rate (m³/d); and i is the stage number.

## Technical bounds

The state-of-the-art recommends keeping the hydraulic loading rate (HLR) between 0.25 m/d (to prevent plant water stress) and 0.7 m/d (to avoid aeration issues). The HLR can be further reduced to 0.1 m/d if a recirculation system is installed (Prost-Boucle & Molle 2012).

Regarding filter layer depths, the minimum 0.3 m, for the filtration layer, and depths exceeding 0.6 m are considered irrelevant (Morvannou et al. 2015; Paing et al. 2015). These values are consistent for both the first and second stages and are detailed in Table S2 in Supplementary Material 3.

All constraints and bounds are summarized in Figure 1 (simplified version; the complete version is available in Supplementary Material 2, Figure S2).

## Optimization

### Inputs and outputs

After defining the model equations, the objective was to optimize the surface area and depth of each stage while meeting the treatment objectives. The resulting system's inputs and outputs are presented in Table 2a and 2b.

### System equations

Optimization involves identifying the optimal solution through an iterative search for the global minimum of an objective function, within predefined constraints. Figure 2 illustrates the system architecture and the interconnection of the key elements.

The objective function is the core equation for optimizing the system. This is the equation for which a minimum must be found, taking into account the constraints and bounds. In optimization problems, the objective function is typically defined so that its minimum corresponds to the optimal solution. This equation was developed to meet two objectives: First, to propose an economically viable solution by minimizing the material volume – and, consequently, the cost of the treatment wetland – while accounting for the cost differences between construction materials (gravel for the first stage and sand for the second stage). Second, to ensure technical viability by maintaining the proposed sizing as close as possible to the optimum CODt loads.

<table><caption>**Table 2** | Inputs (a) and outputs (b) of the methodology</caption><thead><tr><th>a</th><th>Type</th><th>Unit</th><th>b</th><th>Type</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="5">Input</td><td rowspan="4">Inlet concentrations</td><td>[TSS]<sub>in</sub></td><td>mgTSS/L</td><td rowspan="4">Output</td><td rowspan="4">Outlet concentration</td><td>[TSS]<sub>out</sub></td><td>mgTSS/L</td></tr><tr><td>[BOD<sub>5</sub>]<sub>in</sub></td><td>mgBOD<sub>5</sub>/L</td><td>[BOD<sub>5</sub>]<sub>out</sub></td><td>mgBOD<sub>5</sub>/L</td></tr><tr><td>[TKN]<sub>in</sub></td><td>mgTKN/L</td><td>[TKN]<sub>out</sub></td><td>mgTKN/L</td></tr><tr><td>[COD<sub>t</sub>]<sub>in</sub></td><td>mgCOD<sub>t</sub>/L</td><td>[COD<sub>t</sub>]<sub>out</sub></td><td>mgCOD<sub>t</sub>/L</td></tr><tr><td>Hydraulic flow rate</td><td>L/d</td><td></td><td></td><td>Optimized first stage surface area (S<sub>1</sub>)</td><td>m<sup>2</sup></td></tr><tr><td rowspan="4">Objective concentrations</td><td rowspan="2">[TSS]<sub>obj</sub></td><td rowspan="2">mgTSS/L</td><td rowspan="4"></td><td>Optimized second stage surface area (S<sub>2</sub>)</td><td rowspan="4">m<sup>2</sup></td></tr><tr><td>Optimized first stage depth (Z<sub>1</sub>)</td></tr><tr><td>[BOD<sub>5</sub>]<sub>obj</sub></td><td>mgBOD<sub>5</sub>/L</td><td>Optimized second stage depth (Z<sub>2</sub>)</td><td>m</td></tr><tr><td>[TKN]<sub>obj</sub></td><td>mgTKN/L</td><td></td><td>m</td></tr><tr><td></td><td>[COD<sub>t</sub>]<sub>obj</sub></td><td>mgCOD<sub>t</sub>/L</td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Figure 2** | Architecture and key components of our system.</figcaption></figure>

The objective function is:

$$ \argmin_{S_1, Z_1, S_2, Z_2} (\text{Obj}) = S_1 \times Z_1 \times n_1 \times r_{\text{gravel}} + S_2 \times Z_2 \times n_2 \times \frac{r_{\text{sand}}}{r_{\text{gravel}}} + \frac{(L_{\text{opti CODt 1}} - L_{\text{CODt 1}})^2}{L_{\text{opti CODt 1}}^2} \times \alpha \quad (8) $$

where argmin<sub>S<sub>1</sub>,Z<sub>1</sub>,S<sub>2</sub>,Z<sub>2</sub></sub> (Obj) is the objective function value (unitless); S<sub>i</sub> is the surface area of the stage (m<sup>2</sup>); Z<sub>i</sub> is the depth of the stage (m); n<sub>i</sub> is the number of filters in parallel for each stage (unitless); r<sub>gravel</sub> is the cost per cubic meter of gravel (€/m<sup>3</sup>); r<sub>sand</sub> is the cost per cubic meter of sand (€/m<sup>3</sup>); L<sub>opti CODt i</sub> is the total COD surface load for the stage (gCOD<sub>t</sub>/m<sup>2</sup>/d); L<sub>opti CODt i</sub> is the optimal total COD surface load for the stage (gCOD<sub>t</sub>/m<sup>2</sup>/d); α is the alpha hyperparameter (unitless); and i is the stage number.

This equation consists of several components: the outputs to be optimized and the model parameters also referred to as model hyperparameters.

As shown in Table 2b, the outputs represent the values to be minimized while adhering to the system constraints. These values include the surface area (S<sub>1</sub>) and depth (Z<sub>1</sub>) of the first stage, as well as the surface area (S<sub>2</sub>) and depth (Z<sub>2</sub>) of the second stage.

The value L<sub>CODt 1</sub> represents the COD<sub>t</sub> load entering the first filter stage. It is also an output that will be optimized, as it includes the surface area of the first stage. It is essential to maintain this load as close as possible to the optimal load (L<sub>opti CODt 1</sub>) to ensure proper system operation. It is set to the maximum COD<sub>t</sub> load (350 gCOD<sub>t</sub>/m<sup>2</sup>/d, see Table 1) until an output of 12 mgTKN/L is achieved (threshold value, set by expert knowledge). Once the threshold is reached, the optimal load will decrease linearly until it reaches half the optimal load for an output of 6 mgNTK/L (i.e. 175 gCOD<sub>t</sub>/m<sup>2</sup>/d). This approach reduces oxygen demand in the first stage when nitrification objectives are particularly stringent.

The parameters r<sub>gravel</sub> and r<sub>sand</sub> represent the economic cost ratio between gravel and sand. Since sand is generally more expensive, setting a higher r<sub>sand</sub> coefficient than r<sub>gravel</sub> will favor a second stage with a smaller volume compared with the first stage.

α is a model hyperparameter that does not have physical reality. α ensures a fair balance between the value of the COD<sub>t</sub> load and the material cost. The addition of this model hyperparameter enables the control of the weight of the term associated with the COD<sub>t</sub> load. Its proper calibration guides the design toward standard values. The procedure for optimizing its value is described below.

The objective function clearly highlights the four parameters to be optimized: the surface areas and depths of each stage. Subsequently, an optimization algorithm, the covariance matrix adaptation evolution strategy (CMA-ES) (Hansen & Ostermeier 2001), was employed and is presented in the following section. Its goal is to find the optimal values for these four parameters in order to minimize the objective function.
## Covariance matrix adaptation evolution strategy

When addressing a non-convex objective function, it can be challenging to find the best overall solution. Multiple local minima may exist that closely resemble the desired global minimum (Horst et al. 2000). Some optimization algorithms are limited to finding only local minima, while others make assumptions that may not apply to all situations, such as the search space being convex (Nocedal & Wright 2006). To solve our optimization problem, which has no guarantee of being convex, the evolutionary algorithm CMA-ES, first introduced by Hansen & Ostermeier (2001), was employed. CMA-ES is a stochastic optimization algorithm designed for the robust optimization of non-linear and non-convex problems. A key benefit of CMA-ES is its ability to accumulate information across iterations and adjust its parameters, particularly the covariance matrix of the candidate solutions. This adaptability enhances resilience in the presence of noise (Hansen et al. 2010; Hansen 2023). The Pycma library was used to implement CMA-ES (Hansen et al. 2019).

## Model hyperparametrer

Before implementing the optimization, it was essential to determine the optimal value for the $\alpha$ parameter within the objective function. To determine this value, a simultaneous optimization using CMA-ES was performed to obtain a parameter value that satisfies the two sizing rules outlined below. These rules were established using expert knowledge:

* For outlet target concentrations equal to or exceeding 12 mgTKN/L, the sizing should be as close as possible to 1.2 m³/ people equivalent (PE) for the first stage.

* For an outlet target concentration of 6 mgTKN/L, the sizing should approach 1.5 m³/PE for the first stage.

These sizing rules apply to an average influent, with its characteristics detailed in Table 3. Through the strategic setting of inlet and outlet concentrations, we have optimized the coefficient to achieve the desired sizing, aligning with these specific guidelines.

## Examples for benchmarking

It is essential to tailor the various model parameters to align with the specific local context is essential. In this project, we used parameters corresponding to a temperate climate and a mainland France context. Typical values for $r_{gravel}$ and $r_{sand}$ were selected based on local economic conditions. Materials costs (gravel and sand) were sourced from various suppliers. The cost of gravel ranged from €60–€85/m³, while sand costs ranged from €80–€120/m³. The reference values used in this study were €75/m³ for gravel and €100/m³ for sand, representing the corresponding mean values. For the number of filters in parallel, conventional sizing values for temperate climates were used: three cells for the first stage and two cells for the second (Molle et al. 2005). Please refer to Table S3 in Supplementary Material 3 for a comprehensive summary of the values utilized.

Two synthetic datasets were used to benchmark the model: one representing an average influent concentration (Case 1) and the other a heavily concentrated influent (Case 2). Influent concentrations were estimated using a dataset collected by Mercoiret (2009), focusing on agglomerations with fewer than 2,000 PE in mainland France (see Table 3). For Case 1, average concentrations were used, while the concentrations corresponding to the 95% centile influent concentration were used for Case 2.

Flow rates were calculated based on concentrations for each case and a ratio of 60 gDOB₅/d/PE for a 1,000 PE capacity:

$$Q(m^3/d) = 60(gO_2/d) \times 1,000 \times \frac{1}{[BOD_5]_{in}} (L/mgO_2) \quad (9)$$

where Q is the flow rate for 1,000 PE (m³/d), and $[BOD_5]_{in}$ is the inlet concentration (mgBOD₅/L).

<table><caption>**Table 3** | Characteristics of the two tested influents</caption><thead><tr><th></th><th>[TSS]<sub>in</sub></th><th>[BOD<sub>5</sub>]<sub>in</sub></th><th>[TKN]<sub>in</sub></th><th>[COD<sub>t</sub>]<sub>in</sub></th><th>Flow (Q)</th></tr></thead><tbody><tr><td>Unit</td><td>mgTSS/L</td><td>mgBOD<sub>5</sub>/L</td><td>mgTKN/L</td><td>mgCOD<sub>t</sub>/L</td><td>m<sup>3</sup>/d</td></tr><tr><td>Average influent</td><td>288</td><td>265</td><td>67</td><td>646</td><td>226</td></tr><tr><td>95% centile influent</td><td>696</td><td>570</td><td>123</td><td>1,341</td><td>105</td></tr></tbody></table>
The characteristics of these influents are summarized in Table 3.

To better understand how the model responds to desired values, the model outputs were plotted across a range of objective output concentrations (see Table 4). These ranges remain consistent across the two tested cases. As the sizing is based solely on TKN and COD, only these parameters are presented. Objective values are fixed for TSS and BOD$_5$.

## Uncertainty analysis

A Monte Carlo analysis was performed to evaluate the reliability of the optimizations (Metropolis & Ulam 1949; Metropolis et al. 1953). The analysis examined how uncertainties in the empirical data used for model calibration affect the results.

The method is outlined as follows: (i) Randomly vary experimental data values using a Normal distribution to simulate measurement uncertainties; (ii) Perform linear regressions on the new experimental values to obtain updated model equations; and (iii) Integrate the updated model into the optimization system and carry out a new optimization. By repeating this method a thousand times using Monte Carlo simulation, it is possible to obtain a range of uncertainties related to sizing and the various quantiles. This study focuses exclusively on surface area sizing with TKN.

# RESULTS AND DISCUSSION

## Model validation

This section provides a comparison of the model results with those obtained from the INRAE database.

## Total Kjeldahl nitrogen

Comparisons were made between the experimental data from the INRAE database and the corrected models (refitted using a non-linear fitting technique instead of log-transformation, Tedoldi et al. 2025) from Molle et al. (2005, 2008) regarding TKN. Table S1 in Supplementary Material 1 presents the results of the statistical analysis (the coefficients of determination ($R^2$) are 0.984 for the first stage and 0.955 for the second stage). Figure 3(a) illustrates the experimental values and regression models, showing treated load as a function of applied load for the first stage. Figure 3(b) presents the same results for the second stage.

Figure 3(a) and 3(b) compares the results of fitting the power model using the data from Molle et al. (2005, 2008) and the INRAE database. A discrepancy exists between the two fitted models, with the difference being more pronounced for the second stage. The INRAE dataset represents a larger pool of treatment wetlands, operating conditions, and design features and, therefore exhibiting a greater variability. The model partially fails to capture this variability ($R^2$ of the corrected Molle et al. (2005, 2008): 0.4880 (first stage), 0.4970 (second stage); $R^2$ of the power regression with the values from the database: 0.5245 (first stage), 0.7058 (second stage)), as it does not account for all the contributing factors (e.g., the frequency of batches, number of days in operation, temperature, operation quality, and so on). Consequently, we consider our model relevant, especially since the treatment plants where the 24-h composite samplings were taken by Molle et al. (2005, 2008) have a known, reliable design, and the measurements were also conducted under known and reliable conditions. Nevertheless, the discrepancies between the model and the experimental data in our database remain substantial, highlighting the need for a Monte Carlo study to assess the impact of uncertainties.

## Total chemical oxygen demand

The resulting equations from the non-linear fitting of COD$_t$ concentrations versus depth for the first (Equation (10a)) and second (Equation (10b)) treatment stages are provided below:

$$ [\text{COD}_{\text{db}}]_{\text{out } 1} = [\text{COD}_{\text{db}}]_{\text{in } 1} \times e^{-3.136 \times Z_1} \quad (10a) $$

$$ [\text{COD}_{\text{db}}]_{\text{out } 2} = [\text{COD}_{\text{db}}]_{\text{in } 2} \times e^{-7.008 \times Z_2} \quad (10b) $$

<table><caption>**Table 4** | Test intervals and values for objectives concentrations</caption><thead><tr><th></th><th>[TSS]<sub>obj</sub></th><th>[BOD<sub>5</sub>]<sub>obj</sub></th><th>[TKN]<sub>obj</sub></th><th>[COD<sub>t</sub>]<sub>obj</sub></th></tr></thead><tbody><tr><td>Unit</td><td>mgTSS/L</td><td>mgBOD<sub>5</sub>/L</td><td>mgTKN/L</td><td>mgCOD<sub>t</sub>/L</td></tr><tr><td>Value/Interval</td><td>20</td><td>20</td><td>[5; 20]</td><td>[20; 80]</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure 3** | (a and b) Treated load as a function of applied load in gTKN/m²/d for the first stage (a) and the second stage (b) of a French vertical flow wetland. Blue dots: experimental values from the INRAE database (Morvannou et al. 2015); blue curve: obtained by power regression with the values from the database; orange curves: corrected Molle et al. (2008) model for the first stage and corrected Molle et al. (2005) for the second stage.</figcaption></figure>

where $[COD_{db}]_{out,i}$ is the concentration of the dissolved biodegradable COD fraction at the outlet of the stage (mgCOD$_{db}$/L), $[COD_{db}]_{in,i}$ is the concentration of the dissolved biodegradable COD fraction at the inlet of the stage (mgCOD$_{db}$/L), $Z_i$ is the depth of the stage (m), and $i$ is the stage number.

These equations and their use are summarized in Figure 1.

The $R^2$ obtained for the regression models are 0.823 (Equation (10a)) and 0.990 (Equation (10b)). The results of the LOO statistic test give us the following frames for the parameters obtained: 3.136 within [2.498; 3.703] (Equation (10a)) and 7.008 within [6.320; 7.413] (Equation (10b)). This suggests a significant impact from certain data points and underscores the need for further experiments to improve the model's predictive efficacy.

As with TKN, we compared our model with 24-h composite samplings from the INRAE database. The data, clustered by depth, are presented in Supplementary Material 5. However, the INRAE database lacks sufficient variation in depth, as it mainly contains filter layers of 0.6 m for the first stage and 0.3 m for the second stage. Thus, it was not possible to statistically validate the results using the database.

## Benchmark results

The results of the hyperparameter fitting process resulted in a value of $\alpha = 200$.

By using the CMA-ES algorithm, the optimized surface area and depth of each stage can be obtained. The objective now is to create a visual representation of the evolution of these optimized results as a function of the TKN and COD treatment objectives, for surface area (Figure 4) and depth (Figure 5), respectively. The initial focus is on TKN. Figure 4 illustrates the evolution of the sizing parameter linked to TKN treatment, specifically the surface area, for each of the two stages and each of the two tested influents. Figure 4 also demonstrates the evolution of the total volumes, which evolve in line with the surface areas. The depths of each of the two stages are not plotted as they remain constant with the TKN treatment target.

The initial observations on these results are as follows: (i) The higher the TKN treatment objective, the greater the surface areas and volumes required to achieve the desired output. (ii) The total volume is higher for heavily loaded effluents compared with the average one. (iii) For average influents, increasing the surface area of the first stage (in yellow) ensures the desired output at very low concentrations (below 8 mgTKN/L). This is because increasing the surface area of the first stage represents a lower cost in the objective function. For the 95% centile influent, the proposed solution for a very high TKN treatment level involves a compromise between increasing the surface areas of the first (yellow) and second (green) stages to optimize the objective function.
<figure><img src="image_11.png" /><figcaption>**Figure 4** | (a–d) Optimized total treatment wetland surface areas and volume of each of the two stages as a function of the TKN treatment target. (a) and (b) show the surface area optimizations (yellow for the first stage and green for the second stage), for the average influent (a) and the 95% centile influent (b). Dashed lines stand for the sizing recommendation currently in use in mainland France: 0.4 m²/PE per filter in operation, so 1.2 m²/PE for the first stage and 0.8 m²/PE for the second stage (in yellow and green, respectively). (c) and (d) are showing the total volume optimizations, for the average influent (c) and the 95% centile influent (d).</figcaption></figure>

Current sizing in France is based on surface load limits (Table 1) and theoretical pollutant production per PE. This typically results in a design with 0.4 m²/PE per filter in operation (represented by dashed lines in Figure 4(a) and 4(b)). An identical surface area is automatically allocated to the filters of the second stage. For typical treatment objectives (from 20 mgTKN/L to 9 mgNTK/L) and considering three filters in parallel for the first stage, the model proposes a total surface area of 1.3 m²/PE. This exceeds the usual sizing of 1.2 m²/PE for the first stage, which can be explained by the TSS load limit, reached as early as 0.43 m²/PE. The filter in operation cannot be sized below this boundary. However, this is offset by a second stage with a smaller surface area than typically required. Overall, for the average influent, the model proposes a sizing that is not significantly lower than 2 m²/PE (approximately 1.3 m²/PE for the first stage and around 0.65 m²/PE for the second stage). Conversely, from an economic perspective, the proposed sizing is more compelling in scenarios where there is a substantial discrepancy in the volumetric cost of sand and gravel.

As the outlet target concentration of TKN drops below 8.2 mgTKN/L, the surface area required per PE to achieve the desired treatment level increases. At approximately 7 mgTKN/L, the surface area exceeds 0.5 m²/PE per filter in operation. These results demonstrate that current guidelines allow for a high level of nitrification, but would need to be adapted if outlet TKN requirements were below 8 mg/L. Furthermore, these results illustrate that sizing of the second stage could be optimized to reduce construction costs and the system's footprint.

In light of the 95% centile influent, it is clear that an increase in surface area is necessary to accommodate higher loads. The first stage of the process remains constrained by the TSS load until the outlet TKN target is below 8 mg/L. The oversizing of
<figure><img src="image_12.png" /><figcaption>**Figure 5** | (a–d) Optimized total treatment wetland depths and volume of each of the two stages as a function of the CODt treatment target. (a) and (b) are showing the depth optimizations (yellow for the first stage and green for the second stage), for the average influent (a) and the 95% centile influent (b). Dashed lines stand for the sizing recommendation currently in use in mainland France: 0.3 m per filter in operation, both for the first and the second stage (in pink). (c) and (d) show the total volume optimizations, for the average influent (c) and the 95% centile influent (d).</figcaption></figure>

the first stage surface area results in a less TKN-concentrated effluent at the outlet, allowing the second stage to have a slightly lower sizing than the second stage of the average influent. However, for more stringent outlet requirements for TKN, the surface of the second stage must also be increased. The model allows for the quantification of the point at which the system approaches its limits and requires adaptation.

In addition to the TKN analysis, a similar investigation into CODt is presented in Figure 5, illustrating the evolution of the sizing parameter related to CODt treatment, focusing specifically on the depths of each of the two stages for the two tested influents and the corresponding material volumes. The surface areas of each stage are not shown as they remain constant regardless of the CODt treatment target.

The initial observations on these results are as follows: (i) The more restrictive the CODt treatment objective, the greater the depths and volumes required to guarantee the desired output. (ii) The total volume is greater by about 5% for treating the more heavily loaded effluent than the average one. (iii) To reach high levels of treatment, after a phase of increasing depth, the two stages reach a plateau. This is linked to the bounds, which limit the depth of each stage. Consequently, below 35 mgCODt/L for the average influent and below 52 mgCODt/L for the 95% influent percentile, this two-stage treatment train is unable to achieve the desired CODt treatment objective. It will therefore be necessary to select an alternative treatment solution.

With regard to depth evolution, the main observation is that oversizing the filtering layer of materials above 30 cm for each of the filters in parallel (see dashed lines in Figure 5(a) and 5(b)) is unnecessary, except in cases where demanding treatment
levels are sought (below 50 mgCOD/L in the case of conventional influent) or in the presence of a highly loaded influent, like the 95% centile influent.

The approach enhances cost-effective sizing within a standard output scope by ensuring the best distribution of volume at each stage, taking into account material costs. For more stringent output targets, it adjusts as needed, possibly allowing for more intensive solutions.

## Uncertainty analysis results

The Monte Carlo analysis is designed to account for the inherent uncertainties associated with the equations used to represent performance on optimized volumes and dimensions. Figure 6 shows the different results obtained by varying the experimental values using a normal law, organized by quantiles. As previously stated, this study focuses on surface area sizing with TKN.

The results demonstrate that uncertainty on model fitting does not significantly affect the optimized sizing value for discharge levels up to 10 mgTKN/L, for both average and 95% centile influents. As illustrated in the figures, the curves corresponding to four quantiles of all the results of the uncertainty analysis are perfectly superimposed, indicating that there is no variation in the optimization results. For outlet required levels below 10 mgTKN/L, the optimization remains unaffected for the average influent (with only slight variation in the 95th quantile when the treatment target is at its lowest, below 6 mgTKN/L). In contrast, discrepancies between the quantile curves of the uncertainty analysis results are observed for the 95% centile influent. As the level of treatment becomes more stringent, optimization results become more sensitive to model-fitting uncertainties. It can therefore be concluded that the optimization process is not significantly affected by model-fitting uncertainty for normal treatment levels, but less so when aiming for very high treatment levels. This is likely because the models fit particularly well the data for a specific range of input concentrations where most of the observation exists. The uncertainty analysis exhibits that using the model outside the range for which it has been specifically trained may lead to erroneous sizing values. Therefore, using the model for average influents is a more reliable approach. In the future, a more precise model for handling heavy loads will be necessary. Nevertheless, accounting for uncertainties allows the designer to determine when it is preferable to oversize the treatment train design (e.g., by following 75th or 95th quantiles) to ensure specific outlet requirements.

<figure><img src="image_13.png" /><figcaption>**Figure 6** | (a and b) Optimized total volume sizing as a function of the desired outlet TKN concentration. The reference, in hashed lines, corresponds to the sizing with model data, while the various quantiles show the dispersion of the results according to a Monte Carlo analysis. For both average influent (a) and 95% centile influent (b).</figcaption></figure>
An additional Monte Carlo study is provided in Supplementary Material 6, demonstrating that the dispersion of results also increases as inlet concentrations vary. If inlet concentrations are very high, uncertainties in volume sizing also increase, particularly when aiming for very stringent discharge levels.

## CONCLUSION

This work introduces a new approach to sizing French vertical flow wetlands. By defining simple regression models that link TKN to surface sizing and COD to depth sizing, the model employs robust equations that facilitate the independent sizing of each stage of the treatment train. This work introduces a new approach to sizing French vertical flow wetlands. Moreover, by integrating constraint equations and establishing a precise objective function, the entire treatment train can be optimized. This approach allows for meeting discharge criteria while reducing the volume of materials used. It enables more precise design, and, in some cases, reduces construction costs. The optimization process has yielded a design that is tailored to both the characteristics of the influent and current regulatory requirements. This approach differs from traditional methods which rely solely on the treatment plant capacity (population equivalents), as it allows for further refinement by optimizing for specific performance targets. Extensive validation work was carried out on both the models and the optimization, particularly in relation to surface area sizing. These validations confirm the effectiveness and reliability of our approach for systems with influents similar to those in our database. They also provide designers with valuable insight into the optimal design choices. These are connected to constraints on output concentration or, if these are not limiting, to the economic factor.

This method also paves the way for future optimizations and broader applications in wastewater treatment. The framework is versatile and, given the availability of data, additional removal equations can be incorporated to address varying climates and local conditions. The proposed approach can be enhanced by adding other treatment wetland processes, to widen the possibilities in terms of treatment systems and combinations (e.g. recirculation, three-stage or parallel processes, etc.). In addition, this approach can be integrated into a broader decision-support process that incorporates new parameters such as social and environmental concerns.

## ACKNOWLEDGEMENTS

The CARIBSAN project is co-financed by the INTERREG Caribbean program under the European Regional Development Fund, the French Development Agency, and the Martinique and Guadeloupe Water Offices.

## DATA AVAILABILITY STATEMENT

All relevant data are available from an online repository at https://forgemia.inra.fr/reversaal/nature-based-solutions/carib-san/wetlandoptimizer.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

* APHA (American Public Health Association) (2012) *Standard Methods for the Examination of Water and Waste Water*, 22nd edn.
  Washington DC: American Public Health Association, American Water Works Association, and Water Environment Federation.

* Brovelli, A. N., Baechler, S., Rossi, L., Langergraber, G. & Barry, D. A. (2007) 'Coupled flow and hydrogeochemical modelling for design and optimization of horizontal flow constructed wetlands', *Presented at the Second International Symposium on Wetland Pollutant Dynamics and Control (WETPOL)*, Tartu, Estonia, pp. 393–395.

* Defo, C., Kaur, R., Bharadwaj, A., Lal, K. & Kumar, P. (2017) Modelling approaches for simulating wetland pollutant dynamics, *Crit. Rev. Environ. Sci. Technol.*, **47**, 1371–1408. https://doi.org/10.1080/10643389.2017.1350126.

* Dotro, G., Langergraber, G., Molle, P., Nivala, J., Puigagut, J., Stein, O. & Von Sperling, M. (2021) *Treatment Wetlands*. London, UK: IWA Publishing.

* Dubber, D. & Gray, N. F. (2010) Replacement of chemical oxygen demand (COD) with total organic carbon (TOC) for monitoring wastewater treatment performance to minimize disposal of toxic analytical waste, *J. Environ. Sci. Health Part A*, **45**, 1595–1600. https://doi.org/10.1080/10934529.2010.506116.

* García, J., Rousseau, D. P. L., Morató, J., Lesage, E., Matamoros, V. & Bayona, J. M. (2010) Contaminant removal processes in subsurface-flow constructed wetlands: a review, *Crit. Rev. Environ. Sci. Technol.*, **40**, 561–661. https://doi.org/10.1080/10643380802471076.
Gillot, S. & Choubert, J.-M. (2010) Biodegradable organic matter in domestic wastewaters: comparison of selected fractionation techniques, *Water Sci. Technol.*, **62**, 630–639. https://doi.org/10.2166/wst.2010.341.

Giraldi, D., De Michieli Vitturi, M. & Iannelli, R. (2010) FITOVERT: a dynamic numerical model of subsurface vertical flow constructed wetlands, *Environ. Modell. Software*, **25**, 633–640. https://doi.org/10.1016/j.envsoft.2009.05.007.

Grau, P., De Gracia, M., Vanrolleghem, P. A. & Ayesa, E. (2007) A new plant-wide modelling methodology for WWTPs, *Water Res.*, **41**, 4357–4372. https://doi.org/10.1016/j.watres.2007.06.019.

Guillaume-Ruty, S. H. Y., Pueyo-Ros, J., Comas, J. & Forquet, N. (2024) A validation workflow for treatment wetland performance data, *Water Sci. Technol.*, **90**, 598–620. https://doi.org/10.2166/wst.2024.182.

Hansen, N. (2023) The CMA Evolution Strategy: A Tutorial. https://arxiv.org/abs/1604.00772.

Hansen, N. & Ostermeier, A. (2001) Completely derandomized self-adaptation in evolution strategies, *Evol. Comput.*, **9**, 159–195. https://doi.org/10.1162/106365601750190398.

Hansen, N., Auger, A., Ros, R., Finck, S. & Pošík, P. (2010) 'Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009', *Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation*. Presented at the GECCO '10: Genetic and Evolutionary Computation Conference, ACM, Portland Oregon USA, pp. 1689–1696.

Hansen, N., Akimoto, Y. & Baudis, P. (2019) CMA-ES/pycma on Github. Zenodo. DOI:10.5281/zenodo.2559634.

Henze, M., Van Loosdrecht, M. C. M., Ekama, G. A. & Brdjanovic, D. (2008) *Biological Wastewater Treatment: Principles, Modelling and Design*. London, UK: IWA Publishing.

Horst, R., Pardalos, P. M. & Thoai, N. V. (2000) *Introduction to Global Optimization (Nonconvex Optimization and Its Applications*, Vol. 21, 2nd edn. Dordrecht, The Netherlands: Springer.

James, G., Witten, D., Hastie, T. & Tibshirani, R. (2013) *An Introduction to Statistical Learning*, Springer Texts in Statistics. New York, NY: Springer New York.

Kadlec, R. H. (2000) *The inadequacy of first-order treatment wetland models*, *Ecol. Eng.*, **15**, 105–119. https://doi.org/10.1016/S0925-8574(99)00039-7.

Kadlec, R. H. & Knight, R. L. (1996) *Treatment Wetlands*. FL, USA: CRC Press.

Kadlec, R. H. & Wallace, S. D. (2009) *Treatment Wetlands*, 2nd edn. Boca Raton, FL: CRC press.

Langergraber, G. & Šimůnek, J. (2005) Modeling variably saturated water flow and multicomponent reactive transport in constructed wetlands, *Vadose Zone J.*, **4**, 924–938. https://doi.org/10.2136/vzj2004.0166.

Langergraber, G. & Šimůnek, J. (2012) Reactive transport modeling of subsurface flow constructed wetlands using the HYDRUS wetland module, *Vadose Zone J.*, **11**, vzj2011.0104. https://doi.org/10.2136/vzj2011.0104.

Liénard, A. (2010) 'Vertical flow constructed wetlands fed with raw sewage: historical review and recent developments in France', *Presented at the 12th International Conference on Wetland Systems for Water Pollution Control*, Venise, Italy, p. 9.

Meijer, S. C. F. (2004) *Theoretical and Practical Aspects of Modelling Activated Sludge Processes*. Utrecht: Meijer.

Mercoiret, L. (2009) *Domestic Wastewater Characteristics in Rural Areas in France* (in French). Onema, Paris, France.

Metropolis, N. & Ulam, S. (1949) *The Monte Carlo method*, *J. Am. Stat. Assoc.*, **44**, 335–341. https://doi.org/10.1080/01621459.1949.10483310.

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. (1953) Equation of state calculations by fast computing machines, *J. Chem. Phys.*, **21**, 1087–1092. https://doi.org/10.1063/1.1699114.

Meyer, D., Chazarenc, F., Claveau-Mallet, D., Dittmer, U., Forquet, N., Molle, P., Morvannou, A., Pálfy, T., Petitjean, A., Rizzo, A., Samsó Campà, R., Scholz, M., Soric, A. & Langergraber, G. (2015) Modelling constructed wetlands: scopes and aims – a comparative review, *Ecol. Eng.*, **80**, 205–213. https://doi.org/10.1016/j.ecoleng.2014.10.031.

Millot, Y., Troesch, S., Esser, D., Molle, P., Morvannou, A., Gourdon, R. & Rousseau, D. P. L. (2016) Effects of design and operational parameters on ammonium removal by single-stage French vertical flow filters treating raw domestic wastewater, *Ecol. Eng.*, **97**, 516–523. https://doi.org/10.1016/j.ecoleng.2016.10.002.

Ministère de la Transition Écologique (2023) *ROSEAU – Wastewater Treatment Database*. Available at: https://assainissement.developpement-durable.gouv.fr/pages/data/actu.php.

Mitchell, C. & McNevin, D. (2001) Alternative analysis of BOD removal in subsurface flow constructed wetlands employing Monod kinetics, *Water Res.*, **35**, 1295–1303. https://doi.org/10.1016/S0043-1354(00)00373-0.

Molle, P., Liénard, A., Boutin, C., Merlin, G. & Iwema, A. (2005) *How to treat raw sewage with constructed wetlands: an overview of the French systems*, *Water Sci. Technol.*, **51**, 11–21. https://doi.org/10.2166/wst.2005.0277.

Molle, P., Prost-Boucle, S. & Lienard, A. (2008) Potential for total nitrogen removal by combining vertical flow and horizontal flow constructed wetlands: a full-scale experiment study, *Ecol. Eng.*, **34**, 23–29. https://doi.org/10.1016/j.ecoleng.2008.05.016.

Morvannou, A., Choubert, J.-M., Vanclooster, M. & Molle, P. (2014) Modeling nitrogen removal in a vertical flow constructed wetland treating directly domestic wastewater, *Ecol. Eng.*, **70**, 379–386. https://doi.org/10.1016/j.ecoleng.2014.06.034.

Morvannou, A., Forquet, N., Michel, S., Troesch, S. & Molle, P. (2015) Treatment performances of French constructed wetlands: results from a database collected over the last 30 years, *Water Sci. Technol.*, **71**, 1333–1339. https://doi.org/10.2166/wst.2015.089.

Nayak, P. C., Rao, Y. R. S. & Sudheer, K. P. (2006) Groundwater level forecasting in a shallow aquifer using artificial neural network approach, *Water Resour. Manage.*, **20**, 77–90. https://doi.org/10.1007/s11269-006-4007-z.
Naz, M., Uyanik, S., Yesilnacar, M. I. & Sahinkaya, E. (2009) Side-by-side comparison of horizontal subsurface flow and free water surface flow constructed wetlands and artificial neural network (ANN) modelling approach, *Ecol. Eng.*, **35**, 1255–1263. https://doi.org/10.1016/j.ecoleng.2009.05.010.

Nocedal, J. & Wright, S. J. (2006) *Numerical Optimization*, 2nd ed. Springer Series in Operations Research. New York: Springer.

Olsson, L. (2012) *Effect of Design and Dosing Regime on the Treatment Performance of Vertical Flow Constructed Wetlands*. Student thesis, Linköpings universitet.

Paing, J., Guilbert, A., Gagnon, V. & Chazarenc, F. (2015) Effect of climate, wastewater composition, loading rates, system age and design on performances of French vertical flow constructed wetlands: a survey based on 169 full scale systems, *Ecol. Eng.*, **80**, 46–52. https://doi.org/10.1016/j.ecoleng.2014.10.029.

Pálfy, T. G., Molle, P., Langergraber, G., Troesch, S., Gourdon, R. & Meyer, D. (2016) Simulation of constructed wetlands treating combined sewer overflow using HYDRUS/CW2D, *Ecol. Eng.*, **87**, 340–347. https://doi.org/10.1016/j.ecoleng.2015.11.048.

Prost-Boucle, S. & Molle, P. (2012) Recirculation on a single stage of vertical flow constructed wetland: treatment limits and operation modes, *Ecol. Eng.*, **43**, 81–84. https://doi.org/10.1016/j.ecoleng.2012.02.022.

Rieger, L., Takács, I., Villez, K., Siegrist, H., Lessard, P., Vanrolleghem, P. A. & Comeau, Y. (2010) Data reconciliation for wastewater treatment plant simulation studies – planning for high-quality data and typical sources of errors, *Water Environ. Res.*, **82**, 426–433.
https://doi.org/10.2175/106143009X12529484815511.

Rousseau, D. P. L., Vanrolleghem, P. A. & De Pauw, N. (2004) Model-based design of horizontal subsurface flow constructed treatment wetlands: a review, *Water Res.*, **38**, 1484–1493. https://doi.org/10.1016/j.watres.2003.12.013.

Ruiz, H. (2017) *Optimization of French Vertical Flow Treatment Wetlands Applied to Domestic Wastewater Treatment for Different Levels of Performances*, Student thesis, Ecole nationale supérieure Mines-Télécom Atlantique, p. 166.

R Core Team (2024) *R: A language and environment for statistical computing*. R Foundation for Statistical Computing, Vienna, Austria.
https://www.R-project.org/.

Samsó, R. & Garcia, J. (2013a) *BIO_PORE*, a mathematical model to simulate biofilm growth and water quality improvement in porous media: application and calibration for constructed wetlands, *Ecol. Eng.*, **54**, 116–127. https://doi.org/10.1016/j.ecoleng.2013.01.021.

Samsó, R. & García, J. (2013b) Bacteria distribution and dynamics in constructed wetlands based on modelling results, *Sci. Total Environ.*, **461–462**, 430–440. https://doi.org/10.1016/j.scitotenv.2013.04.073.

Stein, O. R., Towler, B. W., Hook, P. B. & Biederman, J. A. (2007) On fitting the k-C\* first order model to batch loaded sub-surface treatment wetlands, *Water Sci. Technol.*, **56**, 93–99. https://doi.org/10.2166/wst.2007.515.

Stone, K. C., Poach, M. E., Hunt, P. G. & Reddy, G. B. (2004) Marsh-pond-marsh constructed wetland design analysis for swine lagoon wastewater treatment, *Ecol. Eng.*, **23**, 127–133. https://doi.org/10.1016/j.ecoleng.2004.07.008.

Sun, G., Cooper, P. F. & Cooper, D. J. (2008) 'The removal of organic matter in horizontal flow reed beds in the UK: performance evaluation using Monod and first-order kinetics', *Presented at the Proceedings of Eleventh International Conference on Wetland Systems for Water Pollution Control*, Indore, India.

Tedoldi, D., Kim, B., Sandoval, S., Forquet, N. & Tassin, B. (2025) Common mistakes and solutions for a better use of correlation- and regression-based approaches in environmental sciences. *Environmental Modelling & Software* **192**, 106526. 10.1016/j.envsoft.2025.106526.

Yuan, C., Huang, T., Zhao, X. & Zhao, Y. (2020) *Numerical models of subsurface flow constructed wetlands: review and future development*, *Sustainability*, **12**, 3498. https://doi.org/10.3390/su12083498.

First received 23 December 2024; accepted in revised form 19 February 2025. Available online 24 May 2025