# What is the best procedure for determining removal rate coefficients in horizontal flow treatment wetlands: influent and effluent concentrations or longitudinal concentration profiles?

Marcos von Sperling, Scott Wallace, Jaime Nivala

► To cite this version:

Marcos von Sperling, Scott Wallace, Jaime Nivala. What is the best procedure for determining removal rate coefficients in horizontal flow treatment wetlands: influent and effluent concentrations or longitudinal concentration profiles?. Water Science and Technology, 2023, 87 (10), pp.2541-2552. 10.2166/wst.2023.144 hal-04220601

HAL Id: hal-04220601

https://hal.science/hal-04220601

Submitted on 18 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
# What is the best procedure for determining removal rate coefficients in horizontal flow treatment wetlands: influent and effluent concentrations or longitudinal concentration profiles?

Marcos von Sperling, Scott Wallace, Jaime Nivala

► To cite this version:

Marcos von Sperling, Scott Wallace, Jaime Nivala. What is the best procedure for determining removal rate coefficients in horizontal flow treatment wetlands: influent and effluent concentrations or longitudinal concentration profiles?. Water Science and Technology, 2023, 87 (10), pp.2541-2552. 10.2166/wst.2023.144 hal-04220601

HAL Id: hal-04220601

https://hal.science/hal-04220601v1

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

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

Water Science & Technology

# What is the best procedure for determining removal rate coefficients in horizontal flow treatment wetlands: influent and effluent concentrations or longitudinal concentration profiles?

Marcos von Sperling IDa, Scott D. Wallaceb and Jaime Nivala Id,c,*

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

ID Mv, 0000-0001-8467-9492; JN, 0000-0003-2990-7688

## ABSTRACT

First-order removal rate coefficients ($k$) are used in predictive equations for estimating effluent concentrations from horizontal flow (HF) wetlands. Due to limited resources, influent and effluent concentration data from existing systems are frequently used in the estimation of $k$ values from operating systems, but another choice is to use concentration data along the longitudinal profile of the HF wetland. Based on a dataset with 41 HF wetlands/studies obtained from a literature survey, with chemical oxygen demand (COD) measurements at different sampling points, volumetric ($k_v$) and areal ($k_A$) removal rate coefficients for the Tanks-In-Series (TIS) model have been obtained using the two estimation methods. In general, removal rate coefficients derived from longitudinal profiles of concentrations were higher than those obtained by using data from influent and effluent concentrations, reflecting the fact that constituent removal is mostly accomplished before the wastewater reaches the outlet zone. Deriving coefficients from longitudinal profiles is more comprehensive, providing a better explanation of the internal removal taking place in the treatment wetland. However, the more widely used approach of calculating $k_v$ and $k_A$ from influent/effluent concentrations may lead to a safer design of horizontal flow wetlands, because of underestimation of the actual removal rate coefficients.

**Key words:** domestic wastewater, first-order kinetics, horizontal subsurface flow constructed wetlands, internal sampling, Tanks-In-Series model

## HIGHLIGHTS

* Lack of information on best way to derive first-order removal rate coefficients.

* Influent/effluent concentrations or longitudinal profile of concentrations are used.

* Database with 41 HF wetlands was constructed and used for calculation of $k$ values.

* Volumetric and areal $k$ values from longitudinal profiles method were usually higher.

* $k$ values from influent/effluent data will result in a more conservative design.

<footnote>aDepartment of Sanitary and Environmental Engineering, Federal University of Minas Gerais (UFMG), Belo Horizonte, Brazil</footnote>

<footnote>bNaturally Wallace Consulting, Pilot Mountain, North Carolina 27041, USA</footnote>

<footnote>cINRAE-French National Research Institute for Agriculture, Food and Environment, Research Unit REVERSAAL, 5 rue de la Doua, CS20244, Villeurbanne Cedex 69625, France</footnote>

<footnote>*Corresponding author. E-mail: jaime.nivala@inrae.fr</footnote>
<figure><img src="image_2.png" /></figure>

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

Water Science & Technology

# What is the best procedure for determining removal rate coefficients in horizontal flow treatment wetlands: influent and effluent concentrations or longitudinal concentration profiles?

Marcos von Sperling IDa, Scott D. Wallaceb and Jaime Nivala Id,c,*

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

ID Mv, 0000-0001-8467-9492; JN, 0000-0003-2990-7688

## ABSTRACT

First-order removal rate coefficients ($k$) are used in predictive equations for estimating effluent concentrations from horizontal flow (HF) wetlands. Due to limited resources, influent and effluent concentration data from existing systems are frequently used in the estimation of $k$ values from operating systems, but another choice is to use concentration data along the longitudinal profile of the HF wetland. Based on a dataset with 41 HF wetlands/studies obtained from a literature survey, with chemical oxygen demand (COD) measurements at different sampling points, volumetric ($k_v$) and areal ($k_A$) removal rate coefficients for the Tanks-In-Series (TIS) model have been obtained using the two estimation methods. In general, removal rate coefficients derived from longitudinal profiles of concentrations were higher than those obtained by using data from influent and effluent concentrations, reflecting the fact that constituent removal is mostly accomplished before the wastewater reaches the outlet zone. Deriving coefficients from longitudinal profiles is more comprehensive, providing a better explanation of the internal removal taking place in the treatment wetland. However, the more widely used approach of calculating $k_v$ and $k_A$ from influent/effluent concentrations may lead to a safer design of horizontal flow wetlands, because of underestimation of the actual removal rate coefficients.

**Key words:** domestic wastewater, first-order kinetics, horizontal subsurface flow constructed wetlands, internal sampling, Tanks-In-Series model

## HIGHLIGHTS

* Lack of information on best way to derive first-order removal rate coefficients.

* Influent/effluent concentrations or longitudinal profile of concentrations are used.

* Database with 41 HF wetlands was constructed and used for calculation of $k$ values.

* Volumetric and areal $k$ values from longitudinal profiles method were usually higher.

* $k$ values from influent/effluent data will result in a more conservative design.

<footnote>aDepartment of Sanitary and Environmental Engineering, Federal University of Minas Gerais (UFMG), Belo Horizonte, Brazil</footnote>

<footnote>bNaturally Wallace Consulting, Pilot Mountain, North Carolina 27041, USA</footnote>

<footnote>cINRAE-French National Research Institute for Agriculture, Food and Environment, Research Unit REVERSAAL, 5 rue de la Doua, CS20244, Villeurbanne Cedex 69625, France</footnote>

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

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

## INTRODUCTION

Removal rate coefficients ($k$) are used in the estimation of the effluent concentrations from biological reactors, in the mathematical description of their performance and also in the comparison of different systems with different configurations or operational modes. In the present analysis, the biological reactor is represented by horizontal subsurface flow treatment wetlands. In most cases for treatment wetlands, first-order reactions are assumed (Kadlec & Wallace 2009), and $k$ is reported as a first-order removal rate coefficient. In saturated horizontal flow (HF) wetlands operating under steady-state conditions, the success in obtaining representative $k$ values depends strongly on a good representation of the reactor hydraulics of the wetland and on a series of practical details that must be considered when running the experiments and processing the data, as listed by von Sperling et al. (2022a). One of these aspects is which sampling points should be used to provide the data on the constituent for which the removal rate coefficient will be calculated: influent and effluent values or values obtained at a profile along the longitudinal axis of the horizontal wetland.

Influent and effluent concentrations are the most widely used data sources in the estimation of $k$ coefficients in treatment wetland studies, regardless of the hydraulic model being used, as specified in a review by Rousseau et al. (2004). This is because these are the most usual sampling or measuring points in any monitoring program, and so existing performance data can be used for the additional benefit of calculating the removal rate coefficient. On the other hand, obtaining longitudinal profile data requires prior planning, since intermediate sampling or measuring points need to be installed at different fractional lengths along the wetland longitudinal axis, and of course, this brings an increase in the number of samples to be analyzed. However, these intermediate points may provide a better representation of the internal behavior of the wetland and, as such, possibly lead to a more representative value of the removal rate coefficient (von Sperling et al. 2020). Figure 1 illustrates these two possibilities, which are the basis for the present study.

As pointed out by von Sperling et al. (2022a), $C_{in}$ and $C_{out}$, or the concentrations measured along the longitudinal profile of the wetland ($C_x$), should be based on average values from a certain monitoring period, and not on single values from one monitoring campaign. This is associated with the influence of the delay related to the liquid traveling time from inlet to outlet, and from one sampling point to the other. To avoid the artifact of transport delay, it is recommended not to compare individual simultaneous measurements at each sampling campaign, but rather averaged values over more than three theoretical hydraulic retention times (Kadlec & Wallace 2009). Long-term removal rate coefficients can be based, for instance, on annual averages or period-of-record (e.g., multiple years of data). As in all sampling campaigns, composite samples are expected to be more accurate than grab samples. Ideally, influent and effluent concentration data should be coupled with inflow and outflow measurements. Average flow values can be calculated using the same averaging strategy employed for the concentrations.

In order to represent mathematically the decay of the constituent concentration along the traveling time or distance from inlet to outlet, a suitable model for the reactor hydraulics needs to be used. The Tanks-In-Series (TIS) model, which is traditionally applied to represent biological reactors in wastewater treatment, is used here. In this hydraulic model, the
## GRAPHICAL ABSTRACT

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

## INTRODUCTION

Removal rate coefficients ($k$) are used in the estimation of the effluent concentrations from biological reactors, in the mathematical description of their performance and also in the comparison of different systems with different configurations or operational modes. In the present analysis, the biological reactor is represented by horizontal subsurface flow treatment wetlands. In most cases for treatment wetlands, first-order reactions are assumed (Kadlec & Wallace 2009), and $k$ is reported as a first-order removal rate coefficient. In saturated horizontal flow (HF) wetlands operating under steady-state conditions, the success in obtaining representative $k$ values depends strongly on a good representation of the reactor hydraulics of the wetland and on a series of practical details that must be considered when running the experiments and processing the data, as listed by von Sperling et al. (2022a). One of these aspects is which sampling points should be used to provide the data on the constituent for which the removal rate coefficient will be calculated: influent and effluent values or values obtained at a profile along the longitudinal axis of the horizontal wetland.

Influent and effluent concentrations are the most widely used data sources in the estimation of $k$ coefficients in treatment wetland studies, regardless of the hydraulic model being used, as specified in a review by Rousseau et al. (2004). This is because these are the most usual sampling or measuring points in any monitoring program, and so existing performance data can be used for the additional benefit of calculating the removal rate coefficient. On the other hand, obtaining longitudinal profile data requires prior planning, since intermediate sampling or measuring points need to be installed at different fractional lengths along the wetland longitudinal axis, and of course, this brings an increase in the number of samples to be analyzed. However, these intermediate points may provide a better representation of the internal behavior of the wetland and, as such, possibly lead to a more representative value of the removal rate coefficient (von Sperling et al. 2020). Figure 1 illustrates these two possibilities, which are the basis for the present study.

As pointed out by von Sperling et al. (2022a), $C_{in}$ and $C_{out}$, or the concentrations measured along the longitudinal profile of the wetland ($C_x$), should be based on average values from a certain monitoring period, and not on single values from one monitoring campaign. This is associated with the influence of the delay related to the liquid traveling time from inlet to outlet, and from one sampling point to the other. To avoid the artifact of transport delay, it is recommended not to compare individual simultaneous measurements at each sampling campaign, but rather averaged values over more than three theoretical hydraulic retention times (Kadlec & Wallace 2009). Long-term removal rate coefficients can be based, for instance, on annual averages or period-of-record (e.g., multiple years of data). As in all sampling campaigns, composite samples are expected to be more accurate than grab samples. Ideally, influent and effluent concentration data should be coupled with inflow and outflow measurements. Average flow values can be calculated using the same averaging strategy employed for the concentrations.

In order to represent mathematically the decay of the constituent concentration along the traveling time or distance from inlet to outlet, a suitable model for the reactor hydraulics needs to be used. The Tanks-In-Series (TIS) model, which is traditionally applied to represent biological reactors in wastewater treatment, is used here. In this hydraulic model, the
<figure><img src="image_6.png" /><figcaption>**Figure 1** | Schematic representation of different sampling strategies for studies of HF wetlands, including the derivation of removal rate coefficients. C is the concentration of the constituent at different positions (influent C<sub>in</sub>, effluent C<sub>out</sub>, and intermediate positions C<sub>x1</sub>, C<sub>x2</sub>, ..., C<sub>xn</sub>). Fractional length (f<sub>x</sub>) is the distance of the sampling point to the inlet (L<sub>i</sub>) divided by the total length (L) of the wetland.</figcaption></figure>

horizontal wetland is conceptually represented by a series of *N* 'apparent' completely mixed tanks. Adaptations of this model have been proposed by Kadlec (2003), who presented the concept of pollutant weathering, by replacing the number of TIS (*N*) with a new parameter (*P*), created to account for both the hydraulics of the reactor and the weathering of a pollutant as it undergoes treatment in the wetland. A further modification was the incorporation of the concept of a background concentration (C*) into the model equations. With *P* and C*, together with the removal rate coefficient *k*, the overall approach was termed P-k-C* (Kadlec & Wallace 2009), which was further widely adopted in the international literature (e.g., textbooks such as Dotro *et al.* 2017; Alarcón Herrera *et al.* 2018; Vidal & Hormazábal 2018).
<figure><img src="image_6.png" /><figcaption>**Figure 1** | Schematic representation of different sampling strategies for studies of HF wetlands, including the derivation of removal rate coefficients. C is the concentration of the constituent at different positions (influent C_in, effluent C_out, and intermediate positions C_x1, C_x2, ..., C xn). Fractional length (f_x) is the distance of the sampling point to the inlet (L_x) divided by the total length (L) of the wetland.</figcaption></figure>

horizontal wetland is conceptually represented by a series of *N* 'apparent' completely mixed tanks. Adaptations of this model have been proposed by Kadlec (2003), who presented the concept of pollutant weathering, by replacing the number of TIS (*N*) with a new parameter (*P*), created to account for both the hydraulics of the reactor and the weathering of a pollutant as it undergoes treatment in the wetland. A further modification was the incorporation of the concept of a background concentration (C*) into the model equations. With *P* and C*, together with the removal rate coefficient *k*, the overall approach was termed P-k-C* (Kadlec & Wallace 2009), which was further widely adopted in the international literature (e.g., textbooks such as Dotro *et al.* 2017; Alarcón Herrera *et al.* 2018; Vidal & Hormazábal 2018).
In spite of the widespread use of the P-k-C* approach, in this paper the traditional structure of the TIS model is used, as recommended by von Sperling et al. (2022b), since it is conceptually more justifiable because it separates reactor hydraulics (as given by N, and not P) from kinetics (as represented by k). The adoption of the values for N and C* are described in the Materials and Methods section of this paper, and the derivation of k is the core of this study.

In the utilization of the TIS models, two approaches can be used:

* Volumetric-based approach: calculations based on the wetland liquid volume (wetland volume × porosity, or V × ε), theoretical hydraulic retention time (HRT or τ), and volumetric removal rate coefficient (k<sub>V</sub>).

* Areal-based approach: calculations based on wetland surface area (A), areal hydraulic loading rate (HLR<sub>A</sub> or q), and areal removal rate coefficient (k<sub>A</sub>).

For design purposes, depending on the approach used (volumetric- or areal-based), the estimation of the effluent concentration (C<sub>out</sub>) from an HF wetland, assuming a first-order reaction, is given by the following equations.

$$C_{\text{out}} = C^* + \frac{(C_{\text{in}} - C^*)}{(1 + k_{\text{V}} \cdot (\tau/N))^N} \quad (1)$$

$$C_{\text{out}} = C^* + \frac{(C_{\text{in}} - C^*)}{(1 + (k_A / N \cdot q))^N} \quad (2)$$

where k<sub>V</sub> is the first-order volumetric removal rate coefficient (d<sup>-1</sup>); k<sub>A</sub> is the first-order areal removal rate coefficient (m/d); C<sub>in</sub> is the influent concentration to wetland (mg/L); C<sub>out</sub> is the effluent concentration from wetland (mg/L); C* is the irreducible background concentration (mg/L); N is the total number of apparent TIS (NTIS) (unitless); τ is the theoretical hydraulic retention time (HRT) of the wetland (d); (τ = V × ε/Q), where V is the volume of the wetland reactor (m<sup>3</sup>), ε is the porosity of the medium (unitless), and Q is the flow rate (m<sup>3</sup>/d); q is the areal hydraulic loading rate, equal to flow divided by surface area [(m<sup>3</sup>/(m<sup>2</sup>·d))].

Areal-based removal rate coefficients (k<sub>A</sub>) and volumetric-based removal rate coefficients (k<sub>V</sub>) can be converted from one another by using the following equation:

$$k_A = k_V \times h \times \varepsilon \quad \text{or} \quad k_V = \frac{k_A}{h \times \varepsilon} \quad (3)$$

where h is the saturated depth (m) and ε is the porosity of the medium (unitless).

The estimation of k<sub>V</sub> coefficients from operating HF wetlands depends on the two methods discussed in this paper: (i) using influent and effluent concentrations (C<sub>in</sub> and C<sub>out</sub>) or (ii) using a longitudinal profile of concentrations along the longitudinal axis of the wetland. The respective equations are given in the following.

## k<sub>V</sub> based on influent and effluent concentrations

The estimation of k<sub>V</sub> coefficients using influent and effluent concentrations is given by Equation (4), which is a rearrangement of Equation (1), making k<sub>V</sub> explicit, and by having measured average values of C<sub>in</sub> and C<sub>out</sub>, assumed values for C*, calculated values of t and estimated values of N.

$$k_V = \frac{N}{\tau} \left[ \left( \frac{C_{\text{in}} - C^*}{C_{\text{out}} - C^*} \right)^{(1/N)} - 1 \right] \quad (4)$$

## k<sub>V</sub> based on the longitudinal profile of concentrations

In this case, there is no explicit equation for making a direct calculation of k<sub>V</sub>. For the longitudinal profile approach, the value of k<sub>V</sub> is obtained by iteration, for instance, using Excel Solver for minimizing the sum of the squared errors related to the observed and estimated concentrations at the different fractional lengths. Therefore, it is necessary to have an equation to represent the decay of the constituent along the longitudinal axis of the HF wetland.

There are two options for calculating the concentrations along the longitudinal profile: smooth profile and stepwise profile. The smooth profile, adopted by Kadlec & Wallace (2009), provides a continuous decaying profile of concentrations along the
In spite of the widespread use of the P-k-C* approach, in this paper the traditional structure of the TIS model is used, as recommended by von Sperling et al. (2022b), since it is conceptually more justifiable because it separates reactor hydraulics (as given by N, and not P) from kinetics (as represented by k). The adoption of the values for N and C* are described in the Materials and Methods section of this paper, and the derivation of k is the core of this study.

In the utilization of the TIS models, two approaches can be used:

* Volumetric-based approach: calculations based on the wetland liquid volume (wetland volume × porosity, or V × ε), theoretical hydraulic retention time (HRT or τ), and volumetric removal rate coefficient (k<sub>V</sub>).

* Areal-based approach: calculations based on wetland surface area (A), areal hydraulic loading rate (HLR<sub>A</sub> or q), and areal removal rate coefficient (k<sub>A</sub>).

For design purposes, depending on the approach used (volumetric- or areal-based), the estimation of the effluent concentration (C<sub>out</sub>) from an HF wetland, assuming a first-order reaction, is given by the following equations.

$$C_{\text{out}} = C^* + \frac{(C_{\text{in}} - C^*)}{(1 + k_{\text{V}} \cdot (\tau / N))^N} \quad (1)$$

$$C_{\text{out}} = C^* + \frac{(C_{\text{in}} - C^*)}{(1 + (k_{\text{A}} / N \cdot q))^N} \quad (2)$$

where k<sub>V</sub> is the first-order volumetric removal rate coefficient (d<sup>-1</sup>); k<sub>A</sub> is the first-order areal removal rate coefficient (m/d); C<sub>in</sub> is the influent concentration to wetland (mg/L); C<sub>out</sub> is the effluent concentration from wetland (mg/L); C* is the irreducible background concentration (mg/L); N is the total number of apparent TIS (NTIS) (unitless); τ is the theoretical hydraulic retention time (HRT) of the wetland (d) (τ = V × ε/Q), where V is the volume of the wetland reactor (m<sup>3</sup>), ε is the porosity of the medium (unitless), and Q is the flow rate (m<sup>3</sup>/d); q is the areal hydraulic loading rate, equal to flow divided by surface area [(m<sup>3</sup>/(m<sup>2</sup>·d))].

Areal-based removal rate coefficients (k<sub>A</sub>) and volumetric-based removal rate coefficients (k<sub>V</sub>) can be converted from one another by using the following equation:

$$k_{\text{A}} = k_{\text{V}} \times h \times \varepsilon \quad \text{or} \quad k_{\text{V}} = \frac{k_{\text{A}}}{h \times \varepsilon} \quad (3)$$

where h is the saturated depth (m) and ε is the porosity of the medium (unitless).

The estimation of k<sub>V</sub> coefficients from operating HF wetlands depends on the two methods discussed in this paper: (i) using influent and effluent concentrations (C<sub>in</sub> and C<sub>out</sub>) or (ii) using a longitudinal profile of concentrations along the longitudinal axis of the wetland. The respective equations are given in the following.

## k<sub>V</sub> based on influent and effluent concentrations

The estimation of k<sub>V</sub> coefficients using influent and effluent concentrations is given by Equation (4), which is a rearrangement of Equation (1), making k<sub>V</sub> explicit, and by having measured average values of C<sub>in</sub> and C<sub>out</sub>, assumed values for C*, calculated values of t and estimated values of N.

$$k_{\text{V}} = \frac{N}{\tau} \left[ \left( \frac{C_{\text{in}} - C^*}{C_{\text{out}} - C^*} \right)^{(1/N)} - 1 \right] \quad (4)$$

## k<sub>V</sub> based on the longitudinal profile of concentrations

In this case, there is no explicit equation for making a direct calculation of k<sub>V</sub>. For the longitudinal profile approach, the value of k<sub>V</sub> is obtained by iteration, for instance, using Excel Solver for minimizing the sum of the squared errors related to the observed and estimated concentrations at the different fractional lengths. Therefore, it is necessary to have an equation to represent the decay of the constituent along the longitudinal axis of the HF wetland.

There are two options for calculating the concentrations along the longitudinal profile: smooth profile and stepwise profile. The smooth profile, adopted by Kadlec & Wallace (2009), provides a continuous decaying profile of concentrations along the
wetland (Equation (5)). The stepwise profile reproduces the conceptual aspect that the wetland is subdivided into N apparent TIS, each tank being completely mixed. Therefore, the stepwise profile of concentrations preserves this concept, and, for a tank with N apparent TIS, there are N different concentrations along the wetland (Equation (6)). Figure 2 illustrates examples of the two estimated profiles. It can be seen that the smooth profile, as the name suggests, leads to a smooth curve, from inlet to outlet. On the other hand, the stepwise profile keeps the same concentration inside each apparent tank in the series, considering the fact that each tank is assumed to be completely mixed. The figure also shows the observed or measured values of the constituent at each fractional length. The values of error in the minimization of the error function for the estimation of $k_V$ are the differences between the observed and estimated values at each fractional length. The minimization algorithm searches for the optimal value of $k_V$, so that the total sum of the squared errors is minimal. In the remainder of this paper, the smooth profile will be used.

Smooth profile:

$$C_x = C^* + \frac{(C_{in} - C^*)}{(1 + k_V(f_x \cdot \tau/N))^N} \quad (5)$$

Stepwise profile:

$$C_{r_n} = C^* + \frac{(C_{in} - C^*)}{(1 + k_V.(\tau/N))^n} \quad (6)$$

where $C_x$ is the estimated concentration at different fractional lengths (mg/L); $f_x$ is the fraction of the longitudinal length along the wetland (e.g., $f_{0.0} = 0.0$ at the inlet, $f_{0.5} = 0.5$ at half of the longitudinal wetland length; $f_{1.0} = 1.0$ at the outlet);

<figure><img src="image_7.png" /><figcaption>**Figure 2** | Example of a plot of longitudinal profile of concentrations in a HF wetland. Circular markers represent the observed or measured concentration values ($C_{observed}$) along the wetland length. Dashed line is the fitted curve for the estimation of concentrations using the smooth profile, according to Equation (5) ($C_{estimated, smooth}$). The solid line is the fitted curve for the estimation of concentrations using the stepwise profile, according to Equation (6) ($C_{estimated, stepwise}$), together with the sequential number of the corresponding apparent tank in the series (numbers below the solid line).</figcaption></figure>
wetland (Equation (5)). The stepwise profile reproduces the conceptual aspect that the wetland is subdivided into N apparent TIS, each tank being completely mixed. Therefore, the stepwise profile of concentrations preserves this concept, and, for a tank with N apparent TIS, there are N different concentrations along the wetland (Equation (6)). Figure 2 illustrates examples of the two estimated profiles. It can be seen that the smooth profile, as the name suggests, leads to a smooth curve, from inlet to outlet. On the other hand, the stepwise profile keeps the same concentration inside each apparent tank in the series, considering the fact that each tank is assumed to be completely mixed. The figure also shows the observed or measured values of the constituent at each fractional length. The values of error in the minimization of the error function for the estimation of $k_V$ are the differences between the observed and estimated values at each fractional length. The minimization algorithm searches for the optimal value of $k_V$, so that the total sum of the squared errors is minimal. In the remainder of this paper, the smooth profile will be used.

Smooth profile:

$$C_x = C^* + \frac{(C_{in} - C^*)}{(1 + k_V(f_x \cdot \tau/N))^N} \quad (5)$$

Stepwise profile:

$$C_{r_n} = C^* + \frac{(C_{in} - C^*)}{(1 + k_V.(\tau/N))^n} \quad (6)$$

where $C_x$ is the estimated concentration at different fractional lengths (mg/L); $f_x$ is the fraction of the longitudinal length along the wetland (e.g., $f_{0.0} = 0.0$ at the inlet, $f_{0.5} = 0.5$ at half of the longitudinal wetland length; $f_{1.0} = 1.0$ at the outlet);

<figure><img src="image_7.png" /><figcaption>**Figure 2** | Example of a plot of longitudinal profile of concentrations in a HF wetland. Circular markers represent the observed or measured concentration values ($C_{observed}$) along the wetland length. Dashed line is the fitted curve for the estimation of concentrations using the smooth profile, according to Equation (5) ($C_{estimated, smooth}$). The solid line is the fitted curve for the estimation of concentrations using the stepwise profile, according to Equation (6) ($C_{estimated, stepwise}$), together with the sequential number of the corresponding apparent tank in the series (numbers below the solid line).</figcaption></figure>
$C_m$ is the concentration at a given apparent reactor along the length in the wetland (mg/L); $r_n$ is the sequential number of the apparent reactor in the series of N reactors (e.g., 1, 2, 3, ..., N) (dimensionless).

Kadlec & Wallace (2009) make an important remark that $k$ values obtained from longitudinal profiles may be higher than those obtained purely from influent/effluent measurements, especially if the effluent concentrations are close to the background concentration $C^*$. The reason is that a concentration close to $C^*$ may be reached before the outlet zone (due to rapid kinetics) and persist close to $C^*$ in the final parts of the wetland. In this case, the actual value of the removal rate coefficient $k$ will be higher than that obtained using measurements of $C_{in}$ and $C_{out}$.

Recognizing the importance of this statement made by Kadlec & Wallace (2009) in terms of deriving $k$ coefficients that may be useful for the design of new HF systems, this paper aims at investigating it in more detail, and with more experimental evidence. Kadlec and Wallace's textbook resorts to a very limited number of references in this particular topic, and most of them are from free water surface wetlands, which have a different hydraulic behavior compared with subsurface flow horizontal wetlands. Therefore, the objective here is to compare both estimation methods (influent/effluent concentrations vs. longitudinal profiles of concentrations) based on experimental evidence from a dataset of several operating HF wetlands in which longitudinal profiles of concentrations were available, together with influent and effluent concentrations.

This comparison of the two estimation methods is not intended to be a mere academic exercise. On a practical level, the relevance of estimating removal rate coefficients in HF wetlands lies in the intended application: (i) generation of $k$ values to be used for the design of future wetlands and (ii) derivation of $k$ values that allow an increase in the understanding of the internal behavior of existing wetlands. These two objectives are discussed in this paper, in syntony with the interpretation of the results obtained with the two estimation methods applied to the dataset of existing wetlands.

## MATERIALS AND METHODS

This study used information available in a dataset originally compiled by Soares (2016) and Ferreira et al. (2020), based on a literature survey of studies that had obtained data on the longitudinal profile of concentrations inside HF wetlands. A total of 41 different wetlands and/or operating conditions have been compiled, and the constituent analyzed was chemical oxygen demand (COD). The dataset, adapted for this study, is presented in the Supplementary Material, and is based on the following references: Valentim (2003), Brasil (2005), Brasil et al. (2007), Sandoval-Cobo & Peña (2007), Chagas (2008), Konnerup et al. (2009), Trang et al. (2010), Chagas et al. (2011), Villaseñor et al. (2011), Eustáquio et al. (2012), and von Sperling & de Paoli (2013). Each study adopted different strategies for sampling frequency and duration of monitoring, and the values of flowrates and influent and effluent concentrations used here for each wetland/study are based on averages from periods of records, as presented by the authors of the publications. From the 41 units/studies investigated, 24.4% received raw influent, 68.3% received primary effluent, and 7.3% received secondary effluent. The descriptive statistics of the main variables are presented in Table 1.

From the table, it can be seen that the HF wetlands used here are mainly at pilot- or small-scale, and the interpretation of the results should take this into account. However, even with this limitation, the dataset is still considered important because, besides providing the internal longitudinal profiles of concentrations, it also includes the physical dimensions of the units (not

<table><caption>**Table 1** | Descriptive statistics of the 41 wetlands/studies that formed the dataset used in this investigation</caption><thead><tr><th>Statistics</th><th>Area A (m²)</th><th>Length L (m)</th><th>Width W (m)</th><th>Liquid depth h (m)</th><th>Aspect ratio (L/W)</th><th>Ratio L/h</th><th>Flow Q (m³/d)</th><th>HRT (d)</th><th>Areal HLR (m³/m²·d)</th><th>Areal OLR (gCOD/m²·d)</th><th>Influent COD C<sub>in</sub> (mg/L)</th><th>Effluent COD C<sub>out</sub> (mg/L)</th><th>Number of tanks in series (N)</th></tr></thead><tbody><tr><td>Minimum</td><td>1.0</td><td>2.0</td><td>0.5</td><td>0.15</td><td>1.0</td><td>3.3</td><td>0.035</td><td>0.8</td><td>0.011</td><td>1.6</td><td>88</td><td>12</td><td>1.5</td></tr><tr><td>Mean</td><td>17.3</td><td>13.6</td><td>1.3</td><td>0.46</td><td>12.3</td><td>53.4</td><td>0.933</td><td>5.3</td><td>0.043</td><td>11.0</td><td>235</td><td>49</td><td>8.5</td></tr><tr><td>Median</td><td>19.2</td><td>12.0</td><td>1.0</td><td>0.55</td><td>7.5</td><td>12.0</td><td>0.448</td><td>4.0</td><td>0.035</td><td>7.7</td><td>197</td><td>33</td><td>3.6</td></tr><tr><td>Stand. dev</td><td>16.2</td><td>10.2</td><td>0.7</td><td>0.25</td><td>10.2</td><td>54.8</td><td>1.466</td><td>5.4</td><td>0.028</td><td>14.7</td><td>130</td><td>40</td><td>7.1</td></tr><tr><td>Maximum</td><td>72.3</td><td>24.1</td><td>4.0</td><td>1.00</td><td>24.0</td><td>160.0</td><td>6.329</td><td>32.0</td><td>0.140</td><td>83.7</td><td>598</td><td>179</td><td>20.7</td></tr></tbody></table>

<footnote>Note: Porosity adopted as 0.35. Number of apparent TIS estimated based on the ratio L/h, according to Equation (7), proposed by von Sperling et al. (2022b): $N = 0.686(L/h)^{0.671}$. Stand. dev., standard deviation; HRT, hydraulic retention time (d); HLR, hydraulic loading rate (q); OLR, organic loading rate.</footnote>
$C_m$ is the concentration at a given apparent reactor along the length in the wetland (mg/L); $r_n$ is the sequential number of the apparent reactor in the series of N reactors (e.g., 1, 2, 3, ..., N) (dimensionless).

Kadlec & Wallace (2009) make an important remark that $k$ values obtained from longitudinal profiles may be higher than those obtained purely from influent/effluent measurements, especially if the effluent concentrations are close to the background concentration $C^*$. The reason is that a concentration close to $C^*$ may be reached before the outlet zone (due to rapid kinetics) and persist close to $C^*$ in the final parts of the wetland. In this case, the actual value of the removal rate coefficient $k$ will be higher than that obtained using measurements of $C_{in}$ and $C_{out}$.

Recognizing the importance of this statement made by Kadlec & Wallace (2009) in terms of deriving $k$ coefficients that may be useful for the design of new HF systems, this paper aims at investigating it in more detail, and with more experimental evidence. Kadlec and Wallace's textbook resorts to a very limited number of references in this particular topic, and most of them are from free water surface wetlands, which have a different hydraulic behavior compared with subsurface flow horizontal wetlands. Therefore, the objective here is to compare both estimation methods (influent/effluent concentrations vs. longitudinal profiles of concentrations) based on experimental evidence from a dataset of several operating HF wetlands in which longitudinal profiles of concentrations were available, together with influent and effluent concentrations.

This comparison of the two estimation methods is not intended to be a mere academic exercise. On a practical level, the relevance of estimating removal rate coefficients in HF wetlands lies in the intended application: (i) generation of $k$ values to be used for the design of future wetlands and (ii) derivation of $k$ values that allow an increase in the understanding of the internal behavior of existing wetlands. These two objectives are discussed in this paper, in syntony with the interpretation of the results obtained with the two estimation methods applied to the dataset of existing wetlands.

## MATERIALS AND METHODS

This study used information available in a dataset originally compiled by Soares (2016) and Ferreira et al. (2020), based on a literature survey of studies that had obtained data on the longitudinal profile of concentrations inside HF wetlands. A total of 41 different wetlands and/or operating conditions have been compiled, and the constituent analyzed was chemical oxygen demand (COD). The dataset, adapted for this study, is presented in the Supplementary Material, and is based on the following references: Valentim (2003), Brasil (2005), Brasil et al. (2007), Sandoval-Cobo & Peña (2007), Chagas (2008), Konnerup et al. (2009), Trang et al. (2010), Chagas et al. (2011), Villaseñor et al. (2011), Eustáquio et al. (2012), and von Sperling & de Paoli (2013). Each study adopted different strategies for sampling frequency and duration of monitoring, and the values of flowrates and influent and effluent concentrations used here for each wetland/study are based on averages from periods of records, as presented by the authors of the publications. From the 41 units/studies investigated, 24.4% received raw influent, 68.3% received primary effluent, and 7.3% received secondary effluent. The descriptive statistics of the main variables are presented in Table 1.

From the table, it can be seen that the HF wetlands used here are mainly at pilot- or small-scale, and the interpretation of the results should take this into account. However, even with this limitation, the dataset is still considered important because, besides providing the internal longitudinal profiles of concentrations, it also includes the physical dimensions of the units (not

<table><caption>**Table 1** | Descriptive statistics of the 41 wetlands/studies that formed the dataset used in this investigation</caption><thead><tr><th>Statistics</th><th>Area A (m²)</th><th>Length L (m)</th><th>Width W (m)</th><th>Liquid depth h (m)</th><th>Aspect ratio (L/W)</th><th>Ratio L/h</th><th>Flow Q (m³/d)</th><th>HRT (d)</th><th>Areal HLR (m³/m²·d)</th><th>Areal OLR (gCOD/m²·d)</th><th>Influent COD C<sub>in</sub> (mg/L)</th><th>Effluent COD C<sub>out</sub> (mg/L)</th><th>Number of tanks in series (N)</th></tr></thead><tbody><tr><td>Minimum</td><td>1.0</td><td>2.0</td><td>0.5</td><td>0.15</td><td>1.0</td><td>3.3</td><td>0.035</td><td>0.8</td><td>0.011</td><td>1.6</td><td>88</td><td>12</td><td>1.5</td></tr><tr><td>Mean</td><td>17.3</td><td>13.6</td><td>1.3</td><td>0.46</td><td>12.3</td><td>53.4</td><td>0.933</td><td>5.3</td><td>0.043</td><td>11.0</td><td>235</td><td>49</td><td>8.5</td></tr><tr><td>Median</td><td>19.2</td><td>12.0</td><td>1.0</td><td>0.55</td><td>7.5</td><td>12.0</td><td>0.448</td><td>4.0</td><td>0.035</td><td>7.7</td><td>197</td><td>33</td><td>3.6</td></tr><tr><td>Stand. dev</td><td>16.2</td><td>10.2</td><td>0.7</td><td>0.25</td><td>10.2</td><td>54.8</td><td>1.466</td><td>5.4</td><td>0.028</td><td>14.7</td><td>130</td><td>40</td><td>7.1</td></tr><tr><td>Maximum</td><td>72.3</td><td>24.1</td><td>4.0</td><td>1.00</td><td>24.0</td><td>160.0</td><td>6.329</td><td>32.0</td><td>0.140</td><td>83.7</td><td>598</td><td>179</td><td>20.7</td></tr></tbody></table>

<footnote>Note: Porosity adopted as 0.35. Number of apparent TIS estimated based on the ratio L/h, according to Equation (7), proposed by von Sperling et al. (2022b): $N = 0.686(L/h)^{0.671}$. Stand. dev., standard deviation; HRT, hydraulic retention time (d); HLR, hydraulic loading rate (q); OLR, organic loading rate.</footnote>
always included in many publications in this field), the hydraulic retention time and applied areal hydraulic loading rate, besides other information.

There was no tracer data with information on the number of apparent TIS (N) in each wetland. Since N is an essential variable for the TIS model and the estimation of the removal rate coefficients ($k_V$ and $k_A$), it was obtained based on empirical Equation (7), proposed by von Sperling et al. (2022b), relating N with the ratio of length/depth. The resulting values are presented in Table 1 and in the Supplementary Material.

$$N = 0.686 \left(\frac{L}{h}\right)^{0.671} \qquad (7)$$

where N is the equivalent number of apparent TIS (dimensionless); L is the wetland length (m); h is the wetland submerged depth (m); and L/h is the length/depth ratio (dimensionless).

For the approach based on influent and effluent concentrations, the volumetric-based removal rate coefficient $k_V$ was estimated using Equation (4) for each wetland/study. For the approach based on longitudinal profile of concentrations, Equation (5) (smooth profile) was used for each wetland/study, together with the minimization of the error function between observed and estimated concentrations at each fractional length using the Solver algorithm within MS Excel. The areal-based removal rate coefficient $k_A$ was calculated according to Equation (3), in which $k_A$ is the product between $k_V$, depth, and porosity.

Porosity was presumed to be 0.35, which is in the range reported in Kadlec & Wallace (2009) for six different studies on HF wetlands (range 0.33–0.41).

In both approaches, $C^*$ was adopted equal to zero. The reason for doing so when calculating removal rate coefficients is that the introduction of $C^*$ in the computations causes an increase in the resulting values of $k_V$ (when applying Equations (4) and (5)). The higher the value of $C^*$, the higher are the resulting $k_V$ and $k_A$ values. Therefore, thinking in terms of design objectives, in order to be on the safe side, it was considered more prudent to adopt a conservative posture and estimate $k_V$ and $k_A$ without considering $C^*$. After that, a designer, when using the predictive equations for estimating effluent concentrations (Equations (1) or (2)), should use the conservative $k_V$ or $k_A$ value, and only then incorporate $C^*$ into the equation for estimating $C_{out}$.

Data on the prevailing liquid temperatures were not available. Therefore, the values of $k_V$ and $k_A$ obtained here are those directly associated with the liquid temperature in each wetland/study, and not the values converted to the standard temperature of 20 °C, as would be, in general, preferable. However, it is believed that this did not impose a limitation with the main objective of the study, which is to compare the estimation methods based on influent/effluent concentrations and longitudinal profile of concentrations, since the liquid temperature in each wetland/study was the same in both calculation approaches.

# RESULTS AND DISCUSSION

## General results of $k_V$ and $k_A$ obtained from both estimation methods

Table 2 presents the resulting volumetric ($k_V$) and areal ($k_A$) removal rate coefficients derived from this study with 41 wetlands/operating conditions, calculated as described above, using the methods based on influent and effluent concentrations and longitudinal profiles of concentrations. The values are expressed as percentiles from the 41 datapoints, in which the 30th and 70th percentiles are meant to represent the lower and upper boundaries of a ‘typical’ range of coefficient values, and the 50th percentile corresponds to the central tendency value.

<table><caption>**Table 2** | Summary of the first-order COD removal rate coefficients derived from the 41 HF wetlands/studies</caption><thead><tr><th rowspan="2">Percentile (%)</th><th colspan="2">Volumetric removal rate coefficient ($k_V$) (d<sup>-1</sup>)</th><th colspan="2">Areal removal rate coefficient ($k_A$) (m/d)</th></tr><tr><th>Based on influent/effluent concentrations</th><th>Based on longitudinal profiles</th><th>Based on influent/effluent concentrations</th><th>Based on longitudinal profiles</th></tr></thead><tbody><tr><td>30</td><td>0.50</td><td>0.58</td><td>0.072</td><td>0.103</td></tr><tr><td><strong>50</strong></td><td><strong>0.68</strong></td><td><strong>0.93</strong></td><td><strong>0.081</strong></td><td><strong>0.125</strong></td></tr><tr><td>70</td><td>0.77</td><td>1.68</td><td>0.104</td><td>0.180</td></tr></tbody></table>

Note: $k_V$ and $k_A$ values are not reported at the standard liquid temperature of 20 °C.
always included in many publications in this field), the hydraulic retention time and applied areal hydraulic loading rate, besides other information.

There was no tracer data with information on the number of apparent TIS (N) in each wetland. Since N is an essential variable for the TIS model and the estimation of the removal rate coefficients ($k_V$ and $k_A$), it was obtained based on empirical Equation (7), proposed by von Sperling et al. (2022b), relating N with the ratio of length/depth. The resulting values are presented in Table 1 and in the Supplementary Material.

$$N = 0.686 \left(\frac{L}{h}\right)^{0.671} \qquad (7)$$

where N is the equivalent number of apparent TIS (dimensionless); L is the wetland length (m); h is the wetland submerged depth (m); and L/h is the length/depth ratio (dimensionless).

For the approach based on influent and effluent concentrations, the volumetric-based removal rate coefficient $k_V$ was estimated using Equation (4) for each wetland/study. For the approach based on longitudinal profile of concentrations, Equation (5) (smooth profile) was used for each wetland/study, together with the minimization of the error function between observed and estimated concentrations at each fractional length using the Solver algorithm within MS Excel. The areal-based removal rate coefficient $k_A$ was calculated according to Equation (3), in which $k_A$ is the product between $k_V$, depth, and porosity.

Porosity was presumed to be 0.35, which is in the range reported in Kadlec & Wallace (2009) for six different studies on HF wetlands (range 0.33–0.41).

In both approaches, $C^*$ was adopted equal to zero. The reason for doing so when calculating removal rate coefficients is that the introduction of $C^*$ in the computations causes an increase in the resulting values of $k_V$ (when applying Equations (4) and (5)). The higher the value of $C^*$, the higher are the resulting $k_V$ and $k_A$ values. Therefore, thinking in terms of design objectives, in order to be on the safe side, it was considered more prudent to adopt a conservative posture and estimate $k_V$ and $k_A$ without considering $C^*$. After that, a designer, when using the predictive equations for estimating effluent concentrations (Equations (1) or (2)), should use the conservative $k_V$ or $k_A$ value, and only then incorporate $C^*$ into the equation for estimating $C_{out}$.

Data on the prevailing liquid temperatures were not available. Therefore, the values of $k_V$ and $k_A$ obtained here are those directly associated with the liquid temperature in each wetland/study, and not the values converted to the standard temperature of 20 °C, as would be, in general, preferable. However, it is believed that this did not impose a limitation with the main objective of the study, which is to compare the estimation methods based on influent/effluent concentrations and longitudinal profile of concentrations, since the liquid temperature in each wetland/study was the same in both calculation approaches.

## RESULTS AND DISCUSSION

### General results of $k_V$ and $k_A$ obtained from both estimation methods

Table 2 presents the resulting volumetric ($k_V$) and areal ($k_A$) removal rate coefficients derived from this study with 41 wetlands/operating conditions, calculated as described above, using the methods based on influent and effluent concentrations and longitudinal profiles of concentrations. The values are expressed as percentiles from the 41 datapoints, in which the 30th and 70th percentiles are meant to represent the lower and upper boundaries of a ‘typical’ range of coefficient values, and the 50th percentile corresponds to the central tendency value.

<table><caption>**Table 2** | Summary of the first-order COD removal rate coefficients derived from the 41 HF wetlands/studies</caption><thead><tr><th rowspan="2">Percentile (%)</th><th colspan="2">Volumetric removal rate coefficient ($k_V$) (d<sup>-1</sup>)</th><th colspan="2">Areal removal rate coefficient ($k_A$) (m/d)</th></tr><tr><th>Based on influent/effluent concentrations</th><th>Based on longitudinal profiles</th><th>Based on influent/effluent concentrations</th><th>Based on longitudinal profiles</th></tr></thead><tbody><tr><td>30</td><td>0.50</td><td>0.58</td><td>0.072</td><td>0.103</td></tr><tr><td><strong>50</strong></td><td><strong>0.68</strong></td><td><strong>0.93</strong></td><td><strong>0.081</strong></td><td><strong>0.125</strong></td></tr><tr><td>70</td><td>0.77</td><td>1.68</td><td>0.104</td><td>0.180</td></tr></tbody></table>

Note: $k_V$ and $k_A$ values are not reported at the standard liquid temperature of 20 °C.
As seen from the table, the removal rate coefficients (both $k_V$ and $k_A$) derived from longitudinal concentration profiles were higher than those that have been obtained using simply influent and effluent concentrations. This topic will be detailed further along in the paper.

Also from Table 2, it is noteworthy to mention that the median areal $k_A$ value for COD (0.081 m/d), obtained from influent/effluent concentrations in this dataset, is comparable to the median value obtained by von Sperling et al. (2022b) for BOD (0.066 m/d, also using influent/effluent data and a completely different dataset). This latter value is based on a separate independent source, comprising an extensive database utilized by Kadlec & Wallace (2009) in their book. The Kadlec and Wallace dataset adopted the P-k-C* approach, which was further converted to the TIS model by von Sperling et al. (2022b).

As mentioned in the Materials and Methods section, the k values reported in Table 2 are not based on the standard liquid temperature of 20 °C, because this information was not available in the sources consulted. This dataset of 41 wetlands/studies is predominantly from warm-climate countries, but the influence of temperature in k values is believed not to be high for the specific case of organic matter, as suggested by Kadlec & Wallace (2009). These authors, working with a large dataset of HF wetlands, using the P-k-C* approach, found the modified-Arrhenius temperature coefficient ($\theta$) for BOD removal to be close to 1.0, varying between a typical range of 0.96 and 1.02 (30th and 70th percentiles). Indeed, Dotro et al. (2017) suggest that for design purposes of new HF wetlands, $\theta$ could be adopted as 1.0 in the case of BOD. This reasoning of minimizing the influence of temperature was also adopted for COD, which is the constituent analyzed in this research.

Regarding the derivation of k values based on the longitudinal profile of concentrations, which are based on the minimization of an error function, the overall fitting between observed and estimated concentrations along the longitudinal distance were very good. The values of the coefficient of determination (CoD) calculated for each of the 41 HF wetlands/studies were systematically high (mean = 0.894; median = 0.916; standard deviation = 0.123 – see Supplementary Material), indicating a good fitting of Equation (5) to represent longitudinal profiles of concentrations (CoD equal to 1.0 indicates a perfect fitting between observed and estimated values). For the influent/effluent method there are no CoD values, since $k_V$ and $k_A$ are directly calculated from the influent and effluent measured concentrations.

## Comparison of $k_V$ and $k_A$ obtained from both estimation methods

Figure 3 (left) presents a comparison of the $k_V$ values obtained using both methods (influent/effluent and longitudinal profile). Looking at the theoretical line with slope 1:1, in which $k_V$ from both methods would be equal, it can be seen that, in general, as already said, $k_V$ based on longitudinal profiles are usually greater than $k_V$ obtained from input and output concentrations. However, it can also be seen that, with this dataset, $k_V$ values lower than $0.6 \text{ d}^{-1}$ were approximately the same with both methods. Figure 3 (right) presents a similar comparison, made for the areal-based $k_A$ coefficient. Also here, with the exception

<figure><img src="image_8.png" /><figcaption>**Figure 3** | Comparison between $k_V$ (left chart) and $k_A$ (right chart) values obtained from the two methods based on influent/effluent COD concentrations and from longitudinal profiles of concentrations, using the dataset of 41 wetlands/studies. The theoretical diagonal line is the 45° line (slope 1:1), where the k values are the same from both methods.</figcaption></figure>

<figure><img src="image_9.png" /></figure>
As seen from the table, the removal rate coefficients (both $k_V$ and $k_A$) derived from longitudinal concentration profiles were higher than those that have been obtained using simply influent and effluent concentrations. This topic will be detailed further along in the paper.

Also from Table 2, it is noteworthy to mention that the median areal $k_A$ value for COD (0.081 m/d), obtained from influent/effluent concentrations in this dataset, is comparable to the median value obtained by von Sperling et al. (2022b) for BOD (0.066 m/d, also using influent/effluent data and a completely different dataset). This latter value is based on a separate independent source, comprising an extensive database utilized by Kadlec & Wallace (2009) in their book. The Kadlec and Wallace dataset adopted the P-k-C* approach, which was further converted to the TIS model by von Sperling et al. (2022b).

As mentioned in the Materials and Methods section, the k values reported in Table 2 are not based on the standard liquid temperature of 20 °C, because this information was not available in the sources consulted. This dataset of 41 wetlands/studies is predominantly from warm-climate countries, but the influence of temperature in k values is believed not to be high for the specific case of organic matter, as suggested by Kadlec & Wallace (2009). These authors, working with a large dataset of HF wetlands, using the P-k-C* approach, found the modified-Arrhenius temperature coefficient (θ) for BOD removal to be close to 1.0, varying between a typical range of 0.96 and 1.02 (30th and 70th percentiles). Indeed, Dotro et al. (2017) suggest that for design purposes of new HF wetlands, θ could be adopted as 1.0 in the case of BOD. This reasoning of minimizing the influence of temperature was also adopted for COD, which is the constituent analyzed in this research.

Regarding the derivation of k values based on the longitudinal profile of concentrations, which are based on the minimization of an error function, the overall fitting between observed and estimated concentrations along the longitudinal distance were very good. The values of the coefficient of determination (CoD) calculated for each of the 41 HF wetlands/studies were systematically high (mean = 0.894; median = 0.916; standard deviation = 0.123 – see Supplementary Material), indicating a good fitting of Equation (5) to represent longitudinal profiles of concentrations (CoD equal to 1.0 indicates a perfect fitting between observed and estimated values). For the influent/effluent method there are no CoD values, since $k_V$ and $k_A$ are directly calculated from the influent and effluent measured concentrations.

## Comparison of $k_V$ and $k_A$ obtained from both estimation methods

Figure 3 (left) presents a comparison of the $k_V$ values obtained using both methods (influent/effluent and longitudinal profile). Looking at the theoretical line with slope 1:1, in which $k_V$ from both methods would be equal, it can be seen that, in general, as already said, $k_V$ based on longitudinal profiles are usually greater than $k_V$ obtained from input and output concentrations. However, it can also be seen that, with this dataset, $k_V$ values lower than $0.6 \text{ d}^{-1}$ were approximately the same with both methods. Figure 3 (right) presents a similar comparison, made for the areal-based $k_A$ coefficient. Also here, with the exception

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

**Figure 3** | Comparison between $k_V$ (left chart) and $k_A$ (right chart) values obtained from the two methods based on influent/effluent COD concentrations and from longitudinal profiles of concentrations, using the dataset of 41 wetlands/studies. The theoretical diagonal line is the 45° line (slope 1:1), where the k values are the same from both methods.

<figure><img src="image_9.png" /></figure>
of the range with low $k_A$ values, in general, the values obtained with the longitudinal profile method are higher than those derived from influent/effluent concentrations. This endorses the statement made by Kadlec & Wallace (2009), the verification of which was one of the main drivers of this study.

In order to understand this a bit further, two actual profiles taken from the 41-wetlands database are shown in Figure 4, aiming at representing two typical situations. The figures plot measured values of COD along the length of two different horizontal wetlands, together with calculated longitudinal profiles using Equation (5) and $k_V$ coefficients obtained employing the two different approaches (measurements of only influent/effluent concentrations and measurements made in different points along the longitudinal axis of the wetland).

The left-hand chart in Figure 4 shows a HF wetland in which the calculation of $k_V$ from both methods led to similar calculated values of COD concentrations and, as such, the two longitudinal profiles are similar, and both have trajectories that are close to the measured values. However, the right-hand graph shows the situation in which $k_V$ is greater when obtained from the longitudinal profile, compared with $k_V$ from influent/effluent concentrations, associated with the fact that the COD removal rate was faster at the initial portions of the wetland. In the right-hand chart, the simulation using $k_V$ from the longitudinal profile shows a better fitting with the measured data, because its calculation was based on minimizing the sum of squared errors (error = observed concentration - estimated concentration) in all measuring points, giving equal weight to each of them. It can be seen that, by searching the best fit for all points has led, in this case, to an under-estimation in the last point (estimated effluent concentration lower than the measured one). For design purposes, it would be a less conservative decision to use this higher $k_V$ value, because it could indicate a lower predicted effluent concentration than was actually obtained in practice. This could be corrected by 'adding in' a background ($C^*$) value when calculating first-order removal rates to be used for design purposes.

However, when analyzing the curve using the lower $k_V$ value, obtained by measured influent/effluent concentrations, it can be seen that the fitting of the curve with the internal points is poor. On the other hand, the effluent concentration is matched exactly – this is expected, because $k_V$ was calculated based on measured $C_{in}$ and $C_{out}$, thus giving guarantee of a perfect fitting to only influent and effluent concentrations. For design purposes, this lower $k_V$ leads to more conservative estimates of effluent concentrations, but underestimates the actual removals within the wetland.

## Applicability of both estimation methods

The sequence of this discussion addresses the points raised in the Introduction section, regarding the applicability of the findings to the design of new HF wetlands, or the investigation of an existing HF wetland. If the wetland designer intends to

<figure><img src="image_10.png" /><figcaption>**Figure 4** | Examples of measured values of COD concentration along the length of two different HF wetlands, together with estimated profiles using $k_V$ coefficients obtained from two different approaches (measurements of only influent/effluent concentrations and measurements made in different points along the longitudinal axis of the wetland). Left graph: $k_V$ values are similar in both methods. Right graph: $k_V$ value calculated from the longitudinal profile is greater than $k_V$ calculated from influent/effluent concentrations.</figcaption></figure>
of the range with low $k_A$ values, in general, the values obtained with the longitudinal profile method are higher than those derived from influent/effluent concentrations. This endorses the statement made by Kadlec & Wallace (2009), the verification of which was one of the main drivers of this study.

In order to understand this a bit further, two actual profiles taken from the 41-wetlands database are shown in Figure 4, aiming at representing two typical situations. The figures plot measured values of COD along the length of two different horizontal wetlands, together with calculated longitudinal profiles using Equation (5) and $k_V$ coefficients obtained employing the two different approaches (measurements of only influent/effluent concentrations and measurements made in different points along the longitudinal axis of the wetland).

The left-hand chart in Figure 4 shows a HF wetland in which the calculation of $k_V$ from both methods led to similar calculated values of COD concentrations and, as such, the two longitudinal profiles are similar, and both have trajectories that are close to the measured values. However, the right-hand graph shows the situation in which $k_V$ is greater when obtained from the longitudinal profile, compared with $k_V$ from influent/effluent concentrations, associated with the fact that the COD removal rate was faster at the initial portions of the wetland. In the right-hand chart, the simulation using $k_V$ from the longitudinal profile shows a better fitting with the measured data, because its calculation was based on minimizing the sum of squared errors (error = observed concentration - estimated concentration) in all measuring points, giving equal weight to each of them. It can be seen that, by searching the best fit for all points has led, in this case, to an under-estimation in the last point (estimated effluent concentration lower than the measured one). For design purposes, it would be a less conservative decision to use this higher $k_V$ value, because it could indicate a lower predicted effluent concentration than was actually obtained in practice. This could be corrected by 'adding in' a background ($C^*$) value when calculating first-order removal rates to be used for design purposes.

However, when analyzing the curve using the lower $k_V$ value, obtained by measured influent/effluent concentrations, it can be seen that the fitting of the curve with the internal points is poor. On the other hand, the effluent concentration is matched exactly – this is expected, because $k_V$ was calculated based on measured $C_{in}$ and $C_{out}$, thus giving guarantee of a perfect fitting to only influent and effluent concentrations. For design purposes, this lower $k_V$ leads to more conservative estimates of effluent concentrations, but underestimates the actual removals within the wetland.

## Applicability of both estimation methods

The sequence of this discussion addresses the points raised in the Introduction section, regarding the applicability of the findings to the design of new HF wetlands, or the investigation of an existing HF wetland. If the wetland designer intends to

<figure><img src="image_10.png" /><figcaption>**Figure 4** | Examples of measured values of COD concentration along the length of two different HF wetlands, together with estimated profiles using $k_V$ coefficients obtained from two different approaches (measurements of only influent/effluent concentrations and measurements made in different points along the longitudinal axis of the wetland). Left graph: $k_V$ values are similar in both methods. Right graph: $k_V$ value calculated from the longitudinal profile is greater than $k_V$ calculated from influent/effluent concentrations.</figcaption></figure>
generate *k* values that will be useful for design purposes, it could be recognized that designers usually err on the side of a more conservative approach. In this case, a conservative approach could be understood as the one that, when there are two estimated values of effluent concentrations from a treatment wetland (inlet and outlet), that method would lead to a higher (in this case, safer) predicted concentration value.

On the other hand, if the objective of estimating the removal rate coefficient *k* is to allow an increase in the understanding of the behavior of an existing wetland, so that the description of this behavior can be useful for researchers and future designs, then the estimation method based on longitudinal profile of concentrations is more useful. This method is associated with a more faithful representation of the internal decrease of the concentration of the constituent under analysis. If the number of studies of concentration profiles increases, allowing the construction of a solid database of *k* values based in internal profiles available in the literature, supported by a large number of HF wetlands with different operating conditions, then this method could start to be the preferred choice. In this case, the designer would decide on how conservative he/she wants to be, based on a better-informed decision procedure.

## Designs for effluent concentrations close to $C^*$

For a designer trying to achieve effluent values close to background concentrations ($C^*$), the two approaches diverge widely. Using the influent/effluent approach, a designer would calculate that a wetland area close to infinity would be required, which is not realistic. An approach using longitudinal profiles with $C^* = 0$ would allow an estimation of the wetland area required to achieve an effluent concentration close to $C^*$ based on the internal profile, which is a more realistic expectation of the technology. It could then be assumed that once the internal profile reaches a value close to $C^*$, further reductions in concentrations would be due to internal biogeochemical cycling within the wetland treatment system. This is a strong benefit of the internal profile approach.

## Lateral homogeneity

A further point deserves additional consideration. In all estimations using profile concentrations, it has been assumed that the values in the longitudinal profile of the HF wetland are a good representation of the whole cross-sectional area in which they have been measured. In reality this may not be the case, due to the existence of preferential flow paths and dead zones, especially in HF wetlands with a large width and poor hydraulic designs in the influent and effluent structures (Wallace & Knight 2006; Fonder et al. 2009). Collecting more samples or making more measurements along the wetland width and then averaging them at each fractional length may be a possibility to reduce this impact, but of course this brings an additional burden in terms of monitoring. However, the difficulties imposed by preferential flow paths and dead zones are not exclusive to longitudinal profiles, since they also affect the estimation of *k* values based on the simpler method of influent and effluent concentrations, because these imperfections will affect the actual mean hydraulic retention time, which will depart from the theoretical hydraulic retention time used in the calculations. A possibility is to consider mean hydraulic retention times given by tracer studies, but this goes beyond the scope of this paper. Nevertheless, it should be considered that the estimated *k* values may incorporate distortions that are not associated to kinetics, but rather to an imperfect hydraulic representation of the reactor.

## Effects of the water balance

Another point that could be brought to this discussion is the influence of the water balance, especially water losses by evapotranspiration, on the interpretation of the measured concentrations. This is addressed by von Sperling et al. (2022a) and is outside the scope of this paper.

## Further discussion on kinetics and hydraulics

These points highlight the difficulties associated with the utilization of simple predictive models for effluent concentrations, such as those with kinetics based on first-order reactions and hydraulics based on an assumed simplified model. Hydraulics and kinetics need to be well characterized in the treatment wetland under study if that data are to be extrapolated to future designs. The TIS model adopted here provides a simple approach to characterizing HF wetland hydraulics, but it is not entirely capable of representing all the intricacies of the internal flows taking place inside the wetland. Furthermore, in terms of kinetics, one would ideally be able to derive and use values of true intrinsic kinetic coefficients. However, this is seldom the case in field-operating treatment wetlands, and the use of 'apparent' removal rate coefficients, still incorporating in themselves the imperfections of the hydraulic model used and other shortcomings, and even the strict assumption of first-
generate *k* values that will be useful for design purposes, it could be recognized that designers usually err on the side of a more conservative approach. In this case, a conservative approach could be understood as the one that, when there are two estimated values of effluent concentrations from a treatment wetland (inlet and outlet), that method would lead to a higher (in this case, safer) predicted concentration value.

On the other hand, if the objective of estimating the removal rate coefficient *k* is to allow an increase in the understanding of the behavior of an existing wetland, so that the description of this behavior can be useful for researchers and future designs, then the estimation method based on longitudinal profile of concentrations is more useful. This method is associated with a more faithful representation of the internal decrease of the concentration of the constituent under analysis. If the number of studies of concentration profiles increases, allowing the construction of a solid database of *k* values based in internal profiles available in the literature, supported by a large number of HF wetlands with different operating conditions, then this method could start to be the preferred choice. In this case, the designer would decide on how conservative he/she wants to be, based on a better-informed decision procedure.

## Designs for effluent concentrations close to $C^*$

For a designer trying to achieve effluent values close to background concentrations ($C^*$), the two approaches diverge widely. Using the influent/effluent approach, a designer would calculate that a wetland area close to infinity would be required, which is not realistic. An approach using longitudinal profiles with $C^* = 0$ would allow an estimation of the wetland area required to achieve an effluent concentration close to $C^*$ based on the internal profile, which is a more realistic expectation of the technology. It could then be assumed that once the internal profile reaches a value close to $C^*$, further reductions in concentrations would be due to internal biogeochemical cycling within the wetland treatment system. This is a strong benefit of the internal profile approach.

## Lateral homogeneity

A further point deserves additional consideration. In all estimations using profile concentrations, it has been assumed that the values in the longitudinal profile of the HF wetland are a good representation of the whole cross-sectional area in which they have been measured. In reality this may not be the case, due to the existence of preferential flow paths and dead zones, especially in HF wetlands with a large width and poor hydraulic designs in the influent and effluent structures (Wallace & Knight 2006; Fonder et al. 2009). Collecting more samples or making more measurements along the wetland width and then averaging them at each fractional length may be a possibility to reduce this impact, but of course this brings an additional burden in terms of monitoring. However, the difficulties imposed by preferential flow paths and dead zones are not exclusive to longitudinal profiles, since they also affect the estimation of *k* values based on the simpler method of influent and effluent concentrations, because these imperfections will affect the actual mean hydraulic retention time, which will depart from the theoretical hydraulic retention time used in the calculations. A possibility is to consider mean hydraulic retention times given by tracer studies, but this goes beyond the scope of this paper. Nevertheless, it should be considered that the estimated *k* values may incorporate distortions that are not associated to kinetics, but rather to an imperfect hydraulic representation of the reactor.

## Effects of the water balance

Another point that could be brought to this discussion is the influence of the water balance, especially water losses by evapotranspiration, on the interpretation of the measured concentrations. This is addressed by von Sperling et al. (2022a) and is outside the scope of this paper.

## Further discussion on kinetics and hydraulics

These points highlight the difficulties associated with the utilization of simple predictive models for effluent concentrations, such as those with kinetics based on first-order reactions and hydraulics based on an assumed simplified model. Hydraulics and kinetics need to be well characterized in the treatment wetland under study if that data are to be extrapolated to future designs. The TIS model adopted here provides a simple approach to characterizing HF wetland hydraulics, but it is not entirely capable of representing all the intricacies of the internal flows taking place inside the wetland. Furthermore, in terms of kinetics, one would ideally be able to derive and use values of true intrinsic kinetic coefficients. However, this is seldom the case in field-operating treatment wetlands, and the use of 'apparent' removal rate coefficients, still incorporating in themselves the imperfections of the hydraulic model used and other shortcomings, and even the strict assumption of first-
order kinetics, is the usual choice that is left to the designer. Nevertheless, the simple approaches adopted here represent a good and effective choice for designers and researchers, provided that they acknowledge the inherent limitations associated with any simple approach to representing the complex environment found in treatment wetlands. Going beyond that would require involving computational fluid dynamics and other categories of advanced mathematical models.

## CONCLUSIONS AND FINAL REMARKS

Based on a database with 41 HF wetlands/studies obtained from the literature, with COD measurements at different sampling points, volumetric ($k_V$) and areal-based ($k_A$) first-order removal rate coefficients for the TIS model have been obtained using two estimation methods associated with different monitoring procedures: (i) influent and effluent concentrations and (ii) longitudinal profiles of concentrations.

It was confirmed that, in general, $k$ coefficients derived from the longitudinal profile of concentrations method tend to be higher than those obtained by using data from only influent and effluent concentrations, indicating that constituent removals may be mostly accomplished before the wastewater reaches the outlet zone.

The median values and typical ranges of areal-based removal rate coefficients ($k_A$) obtained in this study using the influent and effluent concentrations method are comparable to those cited in the literature for organic matter, also based on the same estimation procedure, and representing an independent and larger dataset.

Deriving coefficients from longitudinal profiles is more comprehensive, because it gives a better explanation of the internal removal taking place in the wetland, and so it is expected to be a better descriptor of HF wetlands performance. However, literature is scarce on the availability of $k$ coefficients using longitudinal profiles.

The simple and widely used approach of calculating $k_V$ and $k_A$ from influent/effluent values may lead to a more conservative and safer design of HF wetlands, because the resulting values of the removal rate coefficients are usually lower than those obtained from longitudinal profiles of concentrations. Besides that, most of the $k$ coefficients obtained and reported by researchers are based on this simplified influent/effluent approach.

Although the dataset used here is based mostly on pilot- and small-scale HF wetlands, the results obtained are conceptually coherent and allow an expansion of the understanding of removal rate coefficients in HF wetlands, together with their interpretation, applicability, and limitations.

## ACKNOWLEDGEMENTS

The authors express their gratitude to Dr Alisson Carraro Borges (Federal University of Viçosa, Brazil), for making the original raw dataset of horizontal flow wetlands available. J.N. acknowledges support from the European Commission H2020 Project MULTISOURCE (Grant Agreement 101003527).

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

Alarcón Herrera, M. T., Zurita Martínez, F., Lara-Borrero, J. A. & Vidal, G. 2018 *Humedales de tratamiento: alternativa de saneamiento de aguas residuales aplicable en América Latina*. Pontificia Universidad Javeriana, Bogotá, Colombia, p. 271 (in Spanish).

Brasil, M. S. 2005 *Desempenho de Sistema Alagado Construído para Tratamento de Esgoto Doméstico. [Doutorado em Engenharia Agrícola]*. Departamento de Engenharia Agrícola, Universidade Federal de Viçosa, Viçosa–MG (in Portuguese).

Brasil, M. S., Matos, A. T., Silva, C. M., Cecon, P. R. & Soares, A. A. 2007 Modeling of pollution removal in constructed wetlands with horizontal subsurface flow. *Agric. Eng. Res. 13* (2), 48–56.

Chagas, R. C. 2008 *Utilização de Lírio Amarelo (Hemerocallis flava) em Sistemas Alagados Construídos para tratamento de esgoto doméstico. MSc Dissertation*, Universidade Federal de Viçosa, Viçosa, Minas Gerais, Brasil (in Portuguese).

Chagas, R. C., Matos, A. T., Cecon, P. R., Lo Monaco, P. A. V. & França, L. G. F. 2011 *Cinética de remoção de matéria orgânica em sistemas alagados construídos cultivados com lírio amarelo (Organic matter removal kinetics in constructed wetlands cultivated with yellow lily).* *Rev. Bras. Eng. Agríc. Ambiental-Agriambi* 15 (11), 1186–1192 (in Portuguese).
order kinetics, is the usual choice that is left to the designer. Nevertheless, the simple approaches adopted here represent a good and effective choice for designers and researchers, provided that they acknowledge the inherent limitations associated with any simple approach to representing the complex environment found in treatment wetlands. Going beyond that would require involving computational fluid dynamics and other categories of advanced mathematical models.

## CONCLUSIONS AND FINAL REMARKS

Based on a database with 41 HF wetlands/studies obtained from the literature, with COD measurements at different sampling points, volumetric ($k_V$) and areal-based ($k_A$) first-order removal rate coefficients for the TIS model have been obtained using two estimation methods associated with different monitoring procedures: (i) influent and effluent concentrations and (ii) longitudinal profiles of concentrations.

It was confirmed that, in general, $k$ coefficients derived from the longitudinal profile of concentrations method tend to be higher than those obtained by using data from only influent and effluent concentrations, indicating that constituent removals may be mostly accomplished before the wastewater reaches the outlet zone.

The median values and typical ranges of areal-based removal rate coefficients ($k_A$) obtained in this study using the influent and effluent concentrations method are comparable to those cited in the literature for organic matter, also based on the same estimation procedure, and representing an independent and larger dataset.

Deriving coefficients from longitudinal profiles is more comprehensive, because it gives a better explanation of the internal removal taking place in the wetland, and so it is expected to be a better descriptor of HF wetlands performance. However, literature is scarce on the availability of $k$ coefficients using longitudinal profiles.

The simple and widely used approach of calculating $k_V$ and $k_A$ from influent/effluent values may lead to a more conservative and safer design of HF wetlands, because the resulting values of the removal rate coefficients are usually lower than those obtained from longitudinal profiles of concentrations. Besides that, most of the $k$ coefficients obtained and reported by researchers are based on this simplified influent/effluent approach.

Although the dataset used here is based mostly on pilot- and small-scale HF wetlands, the results obtained are conceptually coherent and allow an expansion of the understanding of removal rate coefficients in HF wetlands, together with their interpretation, applicability, and limitations.

## ACKNOWLEDGEMENTS

The authors express their gratitude to Dr Alisson Carraro Borges (Federal University of Viçosa, Brazil), for making the original raw dataset of horizontal flow wetlands available. J.N. acknowledges support from the European Commission H2020 Project MULTISOURCE (Grant Agreement 101003527).

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

* Alarcón Herrera, M. T., Zurita Martínez, F., Lara-Borrero, J. A. & Vidal, G. 2018 *Humedales de tratamiento: alternativa de saneamiento de aguas residuales aplicable en América Latina*. Pontificia Universidad Javeriana, Bogotá, Colombia, p. 271 (in Spanish).

* Brasil, M. S. 2005 *Desempenho de Sistema Alagado Construído para Tratamento de Esgoto Doméstico. [Doutorado em Engenharia Agrícola].* Departamento de Engenharia Agrícola, Universidade Federal de Viçosa, Viçosa–MG (in Portuguese).

* Brasil, M. S., Matos, A. T., Silva, C. M., Cecon, P. R. & Soares, A. A. 2007 Modeling of pollution removal in constructed wetlands with horizontal subsurface flow. *Agric. Eng. Res. 13* (2), 48–56.

* Chagas, R. C. 2008 *Utilização de Lírio Amarelo (Hemerocallis flava) em Sistemas Alagados Construídos para tratamento de esgoto doméstico. MSc Dissertation*, Universidade Federal de Viçosa, Viçosa, Minas Gerais, Brasil (in Portuguese).

* Chagas, R. C., Matos, A. T., Cecon, P. R., Lo Monaco, P. A. V. & França, L. G. F. 2011 *Cinética de remoção de matéria orgânica em sistemas alagados construídos cultivados com lírio amarelo (Organic matter removal kinetics in constructed wetlands cultivated with yellow lily).* *Rev. Bras. Eng. Agríc. Ambiental-Agriambi* 15 (11), 1186–1192 (in Portuguese).
Dotro, G., Langergraber, G., Molle, P., Nivala, J., Puigagut, J., Stein, O. & Von Sperling, M. 2017 *Treatment Wetlands*, Vol. 7. Biological Wastewater Treatment Series. IWA Publishing, London, UK, p. 154. Available from: https://iwaponline.com/ebooks/book/330/Treatment-Wetlands.

Eustáquio Júnior, V., Matos, A. T., Lo Monaco, P. A. V., Campos, L. C. & Borges, A. C. 2012 Efficiency of constructed wetland systems cultivated with black oats treatment of domestic sewage. *Acta Sci. Technol. Maringá* **34** (4), 391-398.

Ferreira, A. G., Borges, A. C. & Rosa, A. P. 2020 Comparison of first-order kinetic models for sewage treatment in horizontal subsurface-flow constructed wetlands. *Environ. Technol.* doi:10.1080/09593330.2020.1769741.

Fonder, N., Bruch, B., Xanthoulis, D., 2009 Three-dimensional onsite tracer test and the hydraulic behaviour of an operational horizontal subsurface flow treatment wetland. In: *3rd International Symposium on Wetland Pollutant Dynamics and Control, 20-24 September 2009* (Bayona, J. M. & García, J., eds). Universitat Politècnica de Catalunya (UPC), Barcelona, Spain, pp. 177-178.

Kadlec, R. H. 2003 *Effects of pollutant speciation in treatment wetlands design*. Ecol. Eng. **20** (1), 1-16.

Kadlec, R. H. & Wallace, S. D. 2009 *Treatment Wetlands*, 2nd edn. CRC Press, Boca Raton, FL, USA.

Konnerup, D., Koottatep, T. & Brix, H. 2009 *Treatment of domestic wastewater in tropical, subsurface flow constructed wetlands planted with Canna and Heliconia*. Ecol. Eng. **35**, 248-257.

Rousseau, D. P. L., Vanrolleghem, P. A. & Pauw, N. D. 2004 *Model-based design of horizontal subsurface flow constructed treatment wetlands: a review*. Water Res. **38**, 1484-1493.

Sandoval-Cobo, J. J. & Peña, M. R. 2007 Análisis del desempeño de un humedal artificial de flujo sub-superficial en zonas tropicales basado en modelos hidráulicos y una cinética de primer orden. In: *Seminario Manejo Integral de Aguas Residuales Domésticas - Conferencia Latino Americana (LATINOSAN)*, 12-16 November 2007, Cali, Colombia. (in Spanish).

Soares, B. S. 2016 *Comparação de modelos de degradação de matéria orgânica em sistemas alagados construídos tratando esgoto sanitário*. MSc dissertation, Universidade Federal de Viçosa, Brasil (in Portuguese).

Trang, N. T. D., Konnerup, D., Schierup, H. H., Chiem, N. H., Tuan, L. A. & Brix, H. 2010 *Kinetics of pollutant removal from domestic wastewater in a tropical horizontal subsurface flow constructed wetland system: effects of hydraulic loading rate*. Ecol. Eng. **36** (4), 527-535.

Valentim, M. A. A. 2003 *Desempenho de leitos cultivados ('Constructed Wetland') para tratamento de esgoto: contribuições para concepção e operação*. PhD Thesis, Universidade Estadual de Campinas, Campinas, São Paulo, Brasil, pp. 1-233 (in Portuguese).

Vidal, G. & Hormazábal, S. 2018 *Humedales construidos: diseño y operación*. Universidad de Concepción, Concepción, Chile, p. 212 (in Spanish).

Villaseñor, J., Mena, J., Fernández, F. J., Gómez, R. & de Lucas, A. 2011 Kinetics of domestic wastewater COD removal by subsurface flow constructed wetlands using different plant species in temperate period. *Int. J. Environ. Anal. Chem.* **91** (7-8), 693-707.

von Sperling, M. & De Paoli, A. C. 2013 First-order COD decay coefficients associated with different hydraulic models applied to planted and unplanted horizontal subsurface-flow constructed wetlands. *Ecol. Eng.* **57**, 205-209. doi:10.1016/j.ecoleng.2013.04.036.

von Sperling, M., Verbyla, M. E. & Oliveira, S. M. A. C. 2020 *Assessment of Treatment Plant Performance and Water Quality Data: A Guide for Students, Researchers and Practitioners*. IWA Publishing, London, UK. Available at https://iwaponline.com/ebooks/book/777.

von Sperling, M., Wallace, S. D. & Nivala, J. 2022a Practical aspects in the determination of first-order removal rate coefficients (k) for horizontal flow treatment wetlands. In: *Proceedings, IWA 17th International Conference on Wetland Systems for Water Pollution Control*, 6-10 November 2022, Lyon, France.

von Sperling, M., Wallace, S. D. & Nivala, J. 2022b Representing performance of horizontal flow treatment wetlands: the tanks in series (TIS) and the plug flow with dispersion (PFD) approaches and their application to design. *Sci. Total Environ.* **859** (2023), 160259. http://dx.doi.org/10.1016/j.scitotenv.2022.160259.

Wallace, S. & Knight, R. 2006 *Small-Scale Constructed Wetland Treatment Systems; Feasibility, Design Criteria and O&M Requirements*. Water Environment Research Foundation, Alexandria, VA, USA.
Dotro, G., Langergraber, G., Molle, P., Nivala, J., Puigagut, J., Stein, O. & Von Sperling, M. 2017 *Treatment Wetlands*, Vol. 7. Biological Wastewater Treatment Series. IWA Publishing, London, UK, p. 154. Available from: https://iwaponline.com/ebooks/book/330/Treatment-Wetlands.

Eustáquio Júnior, V., Matos, A. T., Lo Monaco, P. A. V., Campos, L. C. & Borges, A. C. 2012 Efficiency of constructed wetland systems cultivated with black oats treatment of domestic sewage. *Acta Sci. Technol. Maringá* **34** (4), 391-398.

Ferreira, A. G., Borges, A. C. & Rosa, A. P. 2020 Comparison of first-order kinetic models for sewage treatment in horizontal subsurface-flow constructed wetlands. *Environ. Technol.* doi:10.1080/09593330.2020.1769741.

Fonder, N., Bruch, B., Xanthoulis, D., 2009 Three-dimensional onsite tracer test and the hydraulic behaviour of an operational horizontal subsurface flow treatment wetland. In: *3rd International Symposium on Wetland Pollutant Dynamics and Control, 20-24 September 2009* (Bayona, J. M. & García, J., eds). Universitat Politècnica de Catalunya (UPC), Barcelona, Spain, pp. 177-178.

Kadlec, R. H. 2003 *Effects of pollutant speciation in treatment wetlands design*. Ecol. Eng. **20** (1), 1-16.

Kadlec, R. H. & Wallace, S. D. 2009 *Treatment Wetlands*, 2nd edn. CRC Press, Boca Raton, FL, USA.

Konnerup, D., Koottatep, T. & Brix, H. 2009 *Treatment of domestic wastewater in tropical, subsurface flow constructed wetlands planted with Canna and Heliconia*. Ecol. Eng. **35**, 248-257.

Rousseau, D. P. L., Vanrolleghem, P. A. & Pauw, N. D. 2004 *Model-based design of horizontal subsurface flow constructed treatment wetlands: a review*. Water Res. **38**, 1484-1493.

Sandoval-Cobo, J. J. & Peña, M. R. 2007 Análisis del desempeño de un humedal artificial de flujo sub-superficial en zonas tropicales basado en modelos hidráulicos y una cinética de primer orden. In: *Seminario Manejo Integral de Aguas Residuales Domésticas - Conferencia Latino Americana (LATINOSAN)*, 12-16 November 2007, Cali, Colombia. (in Spanish).

Soares, B. S. 2016 *Comparação de modelos de degradação de matéria orgânica em sistemas alagados construídos tratando esgoto sanitário*. MSc dissertation, Universidade Federal de Viçosa, Brasil (in Portuguese).

Trang, N. T. D., Konnerup, D., Schierup, H. H., Chiem, N. H., Tuan, L. A. & Brix, H. 2010 *Kinetics of pollutant removal from domestic wastewater in a tropical horizontal subsurface flow constructed wetland system: effects of hydraulic loading rate*. Ecol. Eng. **36** (4), 527-535.

Valentim, M. A. A. 2003 *Desempenho de leitos cultivados ('Constructed Wetland') para tratamento de esgoto: contribuições para concepção e operação*. PhD Thesis, Universidade Estadual de Campinas, Campinas, São Paulo, Brasil, pp. 1-233 (in Portuguese).

Vidal, G. & Hormazábal, S. 2018 *Humedales construidos: diseño y operación*. Universidad de Concepción, Concepción, Chile, p. 212 (in Spanish).

Villaseñor, J., Mena, J., Fernández, F. J., Gómez, R. & de Lucas, A. 2011 Kinetics of domestic wastewater COD removal by subsurface flow constructed wetlands using different plant species in temperate period. *Int. J. Environ. Anal. Chem.* **91** (7-8), 693-707.

von Sperling, M. & De Paoli, A. C. 2013 First-order COD decay coefficients associated with different hydraulic models applied to planted and unplanted horizontal subsurface-flow constructed wetlands. *Ecol. Eng.* **57**, 205-209. doi:10.1016/j.ecoleng.2013.04.036.

von Sperling, M., Verbyla, M. E. & Oliveira, S. M. A. C. 2020 *Assessment of Treatment Plant Performance and Water Quality Data: A Guide for Students, Researchers and Practitioners*. IWA Publishing, London, UK. Available at https://iwaponline.com/ebooks/book/777.

von Sperling, M., Wallace, S. D. & Nivala, J. 2022a Practical aspects in the determination of first-order removal rate coefficients (k) for horizontal flow treatment wetlands. In: *Proceedings, IWA 17th International Conference on Wetland Systems for Water Pollution Control*, 6-10 November 2022, Lyon, France.

von Sperling, M., Wallace, S. D. & Nivala, J. 2022b Representing performance of horizontal flow treatment wetlands: the tanks in series (TIS) and the plug flow with dispersion (PFD) approaches and their application to design. *Sci. Total Environ.* **859** (2023), 160259. http://dx.doi.org/10.1016/j.scitotenv.2022.160259.

Wallace, S. & Knight, R. 2006 *Small-Scale Constructed Wetland Treatment Systems; Feasibility, Design Criteria and O&M Requirements*. Water Environment Research Foundation, Alexandria, VA, USA.