## Abstract

Effluent concentrations from horizontal flow (HF) treatment wetlands can be estimated by using the Tanks-In-Series model for describing hydraulics and first-order removal rate coefficients for describing pollutant removal. In the design of conventional wastewater treatment plants, volumetric removal rate coefficients (*k*_{V}) are traditionally used in conjunction with the theoretical hydraulic retention time. Areal removal rate coefficients (*k*_{A}) coupled with the applied areal hydraulic loading rate are widely used in the literature. Despite this, supporting evidence of its appropriateness is scarce in the literature. The objective of this study is to investigate the adequacy of both approaches by analyzing the influence of liquid depth on *k*_{V} and *k*_{A}. Data from 74 HF wetlands were collected, covering biochemical oxygen demand and chemical oxygen demand, and diverse types of influents (raw sewage and primary, secondary and tertiary effluents). For these conditions, *k*_{V} decreased with depth of the wetland system. Regression analyses between depth and removal rate coefficients were performed, and the equations indicated that *k*_{V} was approximately related to the inverse of depth, while *k*_{A} was almost independent of depth. These findings endorse the utilization of the areal-based approach for design purposes. The volumetric-based approach can also be used, but the value of *k*_{V} must be provided together with the depth being considered.

## HIGHLIGHTS

Literature data were scarce on the adequacy of areal-based approach for horizontal flow wetlands.

Large dataset was used to investigate the influence of depth on removal rate coefficients.

For different types of influents, the volumetric coefficient decreased with wetland depth.

Values of areal removal rate coefficients were found to be little influenced by depth.

Areal-based approach was endorsed by the results of this investigation.

## INTRODUCTION

The present study focuses on horizontal subsurface flow wetlands (HSSF), also characterized by the simpler nomenclature of horizontal flow (HF) wetlands, applied to the treatment of domestic wastewater. The broad topic is the selection of removal rate coefficients to enable the estimation of effluent concentrations at the design stage. In spite of the focus on HF wetlands, this paper uses general concepts of the working principles, kinetics and hydraulics of biological reactors. In the design of conventional wastewater treatment plants, one of the first steps is the estimation of the required dimensions of the biological reactor (length, width, depth). In most cases, the sizing stage starts with the definition of the required reactor volume, based on the criteria available in the general wastewater treatment literature. However, in the case of some treatment processes, and particularly extensive passive systems, such as facultative ponds, overland flow systems, free water surface (FWS) wetlands, HF wetlands and vertical flow wetlands, the surface area – not the volume – is considered the governing factor. The resulting design emphasizes the estimation of the required area (Mara 2003; Kadlec & Wallace 2009; Libhaber & Orozco-Jaramillio 2012; Crites *et al.* 2014). For this, areal hydraulic loading rates and mass loading rates are adopted, based on available local design guidelines, regional experience or specific literature, which, together with information about influent flow and load, lead to the calculation of the required surface area. Reactor depth is then specified, and the resulting volume is calculated based on the area and depth defined at the sizing stage.

After the estimation of the required dimensions for the biological reactor, a subsequent step in the design procedure is the estimation of the expected removal efficiencies and effluent concentrations of the constituents of interest. Frequently, organic matter, expressed in terms of biochemical oxygen demand (carbonaceous BOD_{5}) or chemical oxygen demand (COD), is the selected design target, but nitrogen, phosphorus, coliforms and others may be included, provided that there is sufficient information for the mathematical representation of their removal during treatment. When adopting simple model structures, the estimation of effluent concentrations from most biological reactors in wastewater treatment, including extensive systems, is based on the traditional concept of first-order *volumetric* removal rate coefficients (*k*_{V}) associated with the theoretical *hydraulic retention time* (*τ*) of the reactor. However, in the case of HF treatment wetlands, design and performance evaluation has typically resorted to a different approach, which makes use of the concepts of first-order *areal* removal rate coefficients (*k*_{A}) associated with the *areal hydraulic loading rate* (*q*) applied to the HF wetland. In this adapted way, the surface area of HF wetlands has been recognized as being more influential in treatment performance than the wetland volume itself, not only in terms of sizing of the units but also for the estimation of the effluent concentrations. This adapted focus was the core of most of the design equations presented in the highly influential textbook by Kadlec & Wallace (2009) and has been widely adopted by the international treatment wetlands literature (e.g. textbooks such as Dotro *et al.* 2017; Alarcón Herrera *et al.* 2018; Vidal & Hormazábal 2018).

The equations available in the literature for predicting effluent concentrations according to first-order kinetics and using the traditional volumetric approach incorporate the dimensionless product *k*_{V}·*τ* (*k*_{V} in day^{−1}, *τ* in day), regardless of the hydraulic model utilized for the reactor (completely mixed, plug flow, plug flow with dispersion, tanks in series) (Arceivala 1981; Levenspiel 1999; Metcalf & Eddy 2014). From these equations, it can be seen that the increase in the hydraulic retention time *τ* implies an increase in the removal efficiency of the constituent under analysis. For a certain design inflow rate (*Q*_{in}), the increase in *τ* is obtained by adopting larger reactor volumes (since *τ* = (*V*·*ε*)/*Q*_{in}, where *V* is wetland volume and *ε* is medium porosity), and this can be obtained by increasing the surface area (*A*) and/or the depth (*h*). However, it was reported by Kadlec & Wallace (2009), based on experimental evidence from FWS and HF wetlands, that increasing the volume by increasing the liquid depth was not necessarily accompanied by an enhancement in the pollutant removal rate. On the other hand, when the wetland volume was enlarged by the utilization of greater surface areas in passive wetlands, the positive response in terms of removal rate was more pronounced.

As postulated by Kadlec & Wallace (2009) for these types of passive wetlands, the overall removal rate coefficient (*k*_{V}) may decrease if the depth *h* increases. For instance, as exemplified in their book, for a certain surface area *A*, if *h* is doubled, then *V* doubles, and *τ* is subsequently also doubled. However, the dimensionless product *k*_{V}·*τ* may not double, because the overall *k*_{V} may decrease in a deeper wetland (the deeper zones of a wetland may not be as effective as the shallower zones, which are generally more oxygenated). A possible explanation for HF wetlands is that increasing their liquid level can make the reactor deeper than the rooting depth of the plants, and some functions, such as oxygen diffusion in passive wetlands, are a function of surface area, not depth. Based on this, Kadlec & Wallace (2009) proposed the estimation of removal rates in passive treatment wetlands by the adoption of areal removal rate coefficients (*k*_{A}, expressed in m/day), coupled with the areal hydraulic loading rate (*q*, expressed in m/day or m^{3}/(m^{2}·day), corresponding to inflow *Q*_{in} divided by the wetland surface area *A*). In this case, the dimensionless variable that is an integral part of the first-order removal equations is *k*_{A}/*q*, instead of *k*_{V}·*τ*.

*k*

_{V}) and Equation (2) (based on

*k*

_{A}), for the apparent Tanks-In-Series (TIS) model, which has been widely used in treatment wetlands literature for the estimation of effluent concentrations that are removed according to first-order kinetics (Kadlec & Wallace 2009; Dotro

*et al.*2017):where

*C*

_{in}is the influent concentration to wetland (mg/L);

*C*

_{out}is the effluent concentration from wetland (mg/L);

*C** is the background concentration (mg/L);

*N*is the apparent number of apparent tanks in series (

*N*TIS) (dimensionless);

*k*

_{V}is the volumetric removal rate coefficient (day

^{−1});

*τ*is the total theoretical hydraulic retention time in the wetland (

*τ*=

*V*×

*ε*/

*Q*=

*h*

*×*

*ε*/

*q*);

*Q*= flow (m

^{3}/day);

*V*is the wetland volume (m

^{3});

*h*is the liquid depth (m);

*ε*is the medium porosity (dimensionless);where

*k*

_{A}is the areal removal rate coefficient (m/day);

*q*is the applied areal hydraulic loading rate [m

^{3}/(m

^{2}·day)].

*τ*), liquid volume (

*V*·

*ε*) and flow (

*Q*

_{in}):

*τ*=

*V*·

*ε*/

*Q*

_{in}= (

*A*·

*h)*·

*ε*/

*Q*

_{in}=

*h*·

*ε*/

*q*:

The assumptions that sustain the incisive and important statements that led to the preference of the areal-based approach for removal rate coefficients in the mathematical representation of pollutant removal in HF wetlands obviously require the support of experimental evidence. The book by Kadlec & Wallace (2009) provided, as a main source for HF wetlands, the studies of García *et al.* (2004). Other references on this topic have also been aggregated in this paper. Ideally, this evidence should be based on side-by-side studies, with units in parallel operating with different depths.

García *et al.* (2004), based on pilot-scale HF wetlands operating in parallel (with 0.50 m and with 0.27 m depth, all with the same surface area and receiving the same influent load), analyzed their performance during startup (8 months). All units had the same hydraulic loading rate (*q*), but the shallower ones had approximately half the theoretical hydraulic retention time (*τ*). The results indicated that the beds with the shallower water depth had higher removal efficiencies in terms of COD, BOD, ammonia-N and dissolved reactive phosphorus. In these experiments, *k*_{A} values from the shallower units were higher than those from the deeper units.

Nivala *et al.* (2019) reported results obtained over the course of one year, from two pilot-scale HF wetlands operating in parallel, one with a conventional depth of 0.50 m and the other with 0.25 m. Both had the same surface area, and the deeper wetland received double the flow, so as to keep both units with the same hydraulic retention time. As a result of this strategy, the deeper unit had double the value of the applied area hydraulic loading rate. The effluent concentrations (BOD, total suspended solids, total nitrogen, ammonia-N and *Escherichia coli*) from the shallower unit were predominantly lower than those from the deeper unit, in consonance with the reasoning presented by Kadlec & Wallace (2009). Although not directly included in Nivala *et al.* (2019) paper, a further analysis of their database of results indicated that the *k*_{A} values from both units were similar and, differently from the findings of García *et al.* (2004), the shallower unit did not have higher *k*_{A} values compared with the deeper unit.

Although covering different passive extensive treatment systems, it is worth mentioning that Kadlec & Wallace (2009) describe a study on FWS wetlands, in which the main inference was that increasing depth and hence the hydraulic retention time was not influential in the removal efficiencies of BOD and thermotolerant coliforms. Based on these findings, they postulated the statement that doubling hydraulic retention time by doubling depth reduced the volumetric removal rate *k*_{V} by half, and so the removal efficiency was kept the same, since the dimensionless product *k*_{V}·*τ* ended up being the same. On the other hand, *k*_{A} values remained the same, and this was one of the reasons for advocating the use of *k*_{A} instead of *k*_{V} for passive wetlands. Chiatti & von Sperling (2012) investigated maturation ponds in parallel, with depths of 0.80 and 0.40 m, receiving the same inflow. The results indicated that the shallower ponds, even if having a smaller volume and therefore a shorter retention time, were likely to achieve better nitrogen removal efficiencies than deeper ponds with higher hydraulic retention times. In the case of ponds, the upper layers are determinant in pollutant removal mechanisms, since they represent the highly active photic zone. Based on a similar comparison with maturation ponds in parallel, von Sperling & Mascarenhas (2005) found no statistically significant difference in effluent *E. coli* concentrations from a pond with low depth and low retention time, and another pond in parallel, with double the depth and approximately double the retention time. In this case, disinfection by natural UV radiation, which has a limited penetration in the water column, was the main factor that led to a comparable efficiency in both ponds. Therefore, the additional evidence brought by the analysis of these complementary studies in passive extensive systems support the emphasis than can be placed on areal-based removal mechanisms and rate coefficients.

It can be seen that controlled side-by-side experiments that could cast more light on the comparison between volumetric- and area-based removal rate coefficients are scarce. The objective of this paper is to delve more into the comparison of both removal rate coefficients, but using an alternative method, which makes possible the utilization of data from a large number of operating HF wetlands. The method used here is based on the approach utilized by von Sperling (2005) for facultative and maturation ponds, with a regression analysis between coliform first-order removal rate coefficients with depth, and whose resulting equation clearly indicated that, the lower the pond depth, the higher the first-order volumetric coliform removal rate coefficient. In the present case, data on organic matter removal (BOD and COD) from 74 HF treatment wetlands are used for the regression analysis. The relevance of the study lies mainly on the importance that removal rate coefficients have in the estimation of effluent concentrations from HF treatment wetlands, and that increased understanding about the determination and significance of these coefficients will result in more reliable design calculations for future treatment systems.

## METHODS

In order to obtain a large number of datapoints to perform the regression analyses of depth versus *k*_{V} and *k*_{A}, references from two datasets were used, totaling 74 HF treatment wetlands. The first dataset was originally compiled by Soares (2016) and Ferreira *et al.* (2020), with a total of 41 different HF wetlands and/or operating conditions. The constituent analyzed was 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), von Sperling & de Paoli (2013). From the 41 wetlands/studies investigated, 24% of the systems received raw influent, 60% received primary effluent and 6% received secondary effluent.

The second dataset was obtained from the extensive database used by Kadlec & Wallace (2009) in their analysis of removal rate coefficients from HF treatment wetlands, comprising 182 wetlands/studies. To the author's knowledge, this is the largest dataset on HF wetland performance. Although the database is very large, only 33 wetlands/studies could be used here, because they had information on wetland depth, which is an essential variable for the present study. No data on length and width were available, and the physical characterization of the HF wetlands was based only on surface area and liquid depth. Differently from the first dataset, in this case, the constituent analyzed was BOD. From the 33 units/studies investigated, 46% of the systems received primary effluent, 42% received secondary effluent and 12% received tertiary influent. This dataset of 33 wetlands/studies is also presented in the Supplementary Material.

In order to work with a large database, both datasets have been merged, and Table 1 presents the descriptive statistics of this combined dataset of 74 studies. It can be seen that the dataset presents a large variation in terms of the variables of interest, but with a predominance of small-scale treatment units, as it can be seen from the flows and surface area. Depth, the main input variable for the regression analysis, varied from 0.15 to 1.00 m, although the depth of the systems was largely situated between 0.20 and 0.60 m (10th and 90th percentiles). The median depth was 0.46 m, very close to the value of 0.50 m, considered typical for design purposes, according to Dotro *et al.* (2017). The descriptive statistics of each of the two datasets are presented in the Supplementary Material.

Statistics . | Flow Q (m^{3}/day)
. | Area A (m^{2})
. | Liquid depth h (m)
. | HRT (day) . | Areal HLR q (m^{3}/m^{2}·day)
. |
---|---|---|---|---|---|

Minimum | 0.03 | 1.00 | 0.15 | 0.69 | 0.011 |

10th percentile | 0.10 | 2.00 | 0.20 | 1.48 | 0.022 |

Mean | 4.37 | 80.25 | 0.42 | 4.42 | 0.045 |

Median | 0.66 | 23.25 | 0.46 | 3.75 | 0.035 |

Standard deviation | 17.44 | 299.12 | 0.21 | 4.24 | 0.030 |

90th percentile | 2.29 | 55.00 | 0.60 | 8.82 | 0.088 |

Maximum | 113.60 | 1,750.00 | 1.00 | 32.00 | 0.153 |

Statistics . | Flow Q (m^{3}/day)
. | Area A (m^{2})
. | Liquid depth h (m)
. | HRT (day) . | Areal HLR q (m^{3}/m^{2}·day)
. |
---|---|---|---|---|---|

Minimum | 0.03 | 1.00 | 0.15 | 0.69 | 0.011 |

10th percentile | 0.10 | 2.00 | 0.20 | 1.48 | 0.022 |

Mean | 4.37 | 80.25 | 0.42 | 4.42 | 0.045 |

Median | 0.66 | 23.25 | 0.46 | 3.75 | 0.035 |

Standard deviation | 17.44 | 299.12 | 0.21 | 4.24 | 0.030 |

90th percentile | 2.29 | 55.00 | 0.60 | 8.82 | 0.088 |

Maximum | 113.60 | 1,750.00 | 1.00 | 32.00 | 0.153 |

*Note*: No statistics on organic loading rate and influent and effluent concentrations are presented, because the two databases represented organic matter analysis of two different constituents (BOD and COD). Porosity adopted as 0.35. HRT, theoretical hydraulic retention time; HLR, areal hydraulic loading rate (*q*).

*k*

_{V}was made by Equation (6), which is a rearrangement of Equation (1), making

*k*

_{V}explicit, and by having measured average values of

*C*

_{in}and

*C*

_{out}, assumed values for

*C**, calculated values of

*t*and estimated values of

*N*:

The calculation of the areal removal rate coefficient *k*_{A} was obtained by Equation (4), based on the liquid depth *h* and porosity (). Porosity was adopted as 0.35, which is in the range reported in Kadlec & Wallace (2009) for six different studies on HF wetlands (range 0.33–0.41).

*N*TIS in each wetland to be used in Equation (6). Since

*N*is an essential variable for the estimation of the removal rate coefficients from the TIS model, it was obtained based on empirical Equation (7), proposed by von Sperling

*et al.*(2023), relating

*N*with the ratio of length/depth (

*L/h*). However, this procedure could be employed only for the first dataset, because it contained information on the wetlands’ length. For the second dataset, based on Kadlec & Wallace (2009), no information on length was available, and therefore

*N*could not be estimated for each wetland by this empirical method. The solution was to adopt a fixed value of

*N*for all wetlands/studies in this second dataset, and the value selected was

*N*= 8, which was the median value obtained from tracer studies compiled from a different dataset produced by Kadlec & Wallace (2009) and specified in Dotro

*et al.*(2017) and von Sperling

*et al.*(2023). As will be shown in the Results section, a sensitivity analysis was done, and the impact of the adoption of this fixed value was very small. All

*N*values are presented in the Supplementary Material.where

*N*is the equivalent number of apparent tanks in series (dimensionless);

*L*is the wetland length (m);

*h*is the wetland saturated depth (m);

*L/h*is the length/depth ratio (dimensionless).

The background concentration (*C**) was adopted equal to zero in Equation (6), as suggested by von Sperling *et al.* (2022) in order to lead to more conservative *k* values, which can be more useful for design purposes. Besides that, the first dataset had COD data, while the second dataset worked with BOD concentrations. This was an additional reason that, in order to standardize, all background concentrations were adopted equal to zero.

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. It is expected that this limitation has not brought important implications, based on the suggestions made by Dotro *et al.* (2017), that the modified Arrhenius temperature coefficient (*q*) can be assumed, in the case or organic matter in the design of HF treatment wetlands, equal to 1.0. It can be seen that structuring such a database suffers from many difficulties, related to the diverse nature in which each study in the literature has been conducted. In order to arrive at a common base and overcome the existing limitations, simplifications needed to be introduced. All effort was done to minimize their impact, as explained above, and as will be seen in the Results section.

## RESULTS

### Descriptive statistics of volumetric (*k*_{V}) and areal (*k*_{A}) removal rate coefficients

*k*

_{V}and

*k*

_{A}values obtained for each type of influent (raw sewage, effluent from primary treatment, effluent from secondary treatment, effluent from tertiary treatment). The descriptive statistics of

*k*

_{V}and

*k*

_{A}for each of the two datasets (COD and BOD) are presented separately in the Supplementary Material. Table 2 presents the descriptive statistics (measures of central tendency and dispersion) of

*k*

_{V}and

*k*

_{A}related to the overall dataset with 74 HF treatment wetlands, and Figure 1 complements this information, with percentile values in the form of a box plot. Although BOD and COD are grouped together here, it can be seen that the general results are coherent:

*k*

_{V}and

*k*

_{A}values are higher for raw sewage, which contains a greater fraction of biodegradable organic matter, and are lower for tertiary influents, which are expected to be composed by a larger refractory fraction. Primary and secondary influents are similar and are placed in-between the interval of ‘raw’ and ‘tertiary’. The resulting median areal

*k*

_{A}coefficient for primary and secondary influents is 0.066 m/day, and, interestingly, has exactly the same median value (0.066 m/day) derived from 103 HF wetlands fed by primary and secondary influents, available in the 182 wetlands/studies from the extensive Kadlec & Wallace (2009) database, which was adapted by von Sperling

*et al.*(2023) for the TIS model. Additionally, the range considered typical for design purposes, comprised by the 30th and 70th percentiles for

*k*

_{A}found here (not shown in Table 2), is 0.047 and 0.081 m/day, which is very similar to the one calculated by von Sperling

*et al.*(2022), using the 103 datapoints for primary and secondary influents from Kadlec & Wallace (2009) database (0.048 and 0.100 m/day). From this 103 Kadlec and Wallace datapoints, only 29 are used here, which indicates a satisfactory degree of diversity from both databases and emphasizes the good convergence to similar values of removal rate coefficients.

Statistics . | Volumetric removal rate coefficient k_{V} (day^{−1}). | Areal removal rate coefficient k_{A} (m/day). | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Raw . | Primary . | Secondary . | Tertiary . | Overall . | Raw . | Primary . | Secondary . | Tertiary . | Overall . | |

Number of data (n) | 10 | 43 | 17 | 4 | 74 | 10 | 43 | 17 | 4 | 74 |

Minimum | 0.54 | 0.06 | 0.32 | 0.22 | 0.06 | 0.054 | 0.021 | 0.034 | 0.035 | 0.021 |

Mean | 0.93 | 0.53 | 0.60 | 0.34 | 0.59 | 0.115 | 0.065 | 0.065 | 0.045 | 0.071 |

Median | 0.79 | 0.44 | 0.66 | 0.38 | 0.54 | 0.116 | 0.066 | 0.065 | 0.040 | 0.066 |

Standard deviation | 0.37 | 0.36 | 0.22 | 0.08 | 0.35 | 0.042 | 0.034 | 0.020 | 0.013 | 0.036 |

Maximum | 1.68 | 1.34 | 0.93 | 0.40 | 1.68 | 0.177 | 0.142 | 0.119 | 0.064 | 0.177 |

Statistics . | Volumetric removal rate coefficient k_{V} (day^{−1}). | Areal removal rate coefficient k_{A} (m/day). | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Raw . | Primary . | Secondary . | Tertiary . | Overall . | Raw . | Primary . | Secondary . | Tertiary . | Overall . | |

Number of data (n) | 10 | 43 | 17 | 4 | 74 | 10 | 43 | 17 | 4 | 74 |

Minimum | 0.54 | 0.06 | 0.32 | 0.22 | 0.06 | 0.054 | 0.021 | 0.034 | 0.035 | 0.021 |

Mean | 0.93 | 0.53 | 0.60 | 0.34 | 0.59 | 0.115 | 0.065 | 0.065 | 0.045 | 0.071 |

Median | 0.79 | 0.44 | 0.66 | 0.38 | 0.54 | 0.116 | 0.066 | 0.065 | 0.040 | 0.066 |

Standard deviation | 0.37 | 0.36 | 0.22 | 0.08 | 0.35 | 0.042 | 0.034 | 0.020 | 0.013 | 0.036 |

Maximum | 1.68 | 1.34 | 0.93 | 0.40 | 1.68 | 0.177 | 0.142 | 0.119 | 0.064 | 0.177 |

### Influence of depth on the volumetric removal rate coefficient (*k*_{V}), separated by dataset (BOD and COD) and influent type

*k*

_{V}) as a function of depth. The graphs are separated by dataset, because the first one (41 datapoints, left-hand chart) represented organic matter by COD, and the second one (33 datapoints, right-hand chart) was based on BOD. In each graph, the

*k*

_{V}values are separated by the type of influent to the HF wetlands (raw, primary, secondary and tertiary). In addition, trend lines (power function) were added, to assist in the visualization of the influence of depth on

*k*

_{V}. No regression equations are presented at this stage, because it is considered that there are not enough datapoints to sustain separate regression equations for each influent type. The purpose here is only to analyze the main trends. A regression equation is proposed in the next section, but for the overall dataset (74 datapoints together).

The data from 41 wetlands/studies (left chart) had higher *k*_{V} values, even though they were for COD, because the influents were comprised in large part by raw sewage and primary effluents, which have a higher organic matter biodegradability. The 33 wetlands/studies (right chart) had lower *k*_{V} values, even though they were for BOD, because the influents were of primary, secondary or tertiary origin.

However, what is considered remarkable in Figure 2 is the consistency that all trendlines, for COD and BOD and for all influent types, despite the scatter of data, indicate a decrease of *k*_{V} with depth, supporting the main driver of investigation in this study and the general statement formulated by Kadlec & Wallace (2009).

As pointed out in the Introduction section, it is understandable that the value of *k*_{V} decreases as *h* increases: the wetland zone with higher activity is close to the surface, in the zone of high rhizosphere density. Increasing the depth of the wetland will not improve the beneficial participation of the microbially active upper zone (e.g. rhizosphere). While there are presumably lower levels of microbial activity in the bottom centimeters of HF wetlands, the hydraulic conductivity is likely to be higher, due to a lower density of plant roots and microbes. Shallower HF wetlands (e.g. total saturated depth of 20–30 cm) will force all flow to pass through the root zone, closer to where microbial and rhizospheric activities are highest. Based on microbial community studies at mesocosm scale at different wetland depths (0.10, 0.30 and 0.60 m), Weber & Legge (2013) stated that the catabolic capabilities of the rhizospheric microbial communities were much greater than those of the gravel-associated communities, and that a decrease in catabolic capability was observed with increasing depth, suggesting that microbial communities at shallower depths play a larger role in the degradation of organic matter. García *et al.* (2004), in their experimental studies with different depth HF wetlands, mentioned that the higher efficiency observed in the shallower beds was related to their less reducing conditions, measured by the average redox potential, compared with the deeper beds, what seemed to lead to differences in the biochemical processes. This is similar to the findings of Ren *et al.* (2011), comparing HF wetlands of 0.10, 0.30 and 0.60 m depth, who observed that the highest concentration of dissolved oxygen and the highest amount of plant root production were achieved in the shallower units.

### Regression of depth versus *k*_{V} and *k*_{A} using the overall dataset

*h*) with

*k*

_{V}and

*k*

_{A}, all 74 datapoints were put together (BOD and COD), leading to the regression analyses presented in Figure 3. In this figure, regression equations are presented, because the number of data was considered sufficient. The regression model selected was based on the power function, with the form

*y*=

*a*·

*x*

^{−b}, where

*y*is the dependent variable,

*x*is the independent variable and ‘

*a*’ and ‘

*b*’ are regression coefficients. Other functions were tried for the regression analysis, but the power function showed to give a good fit to the datapoints and also allowed for simple conversions between

*k*

_{V}and

*k*

_{A}, as discussed in the next section.

*h*and

*k*

_{V}(left chart in Figure 3), a wide spread of datapoints can be seen, what could be already expected, since the scatter plot is comprised of wetlands/studies representing different constituents (BOD and COD), influent types (raw, primary, secondary, tertiary), applied loading rates and operating conditions. This is reflected in the not particularly high coefficient of determination (CoD) (0.4204). However, the trend of

*k*

_{V}decreasing with the increase in the depth is clearly seen.

On the other hand, in the regression analysis between *h* and *k*_{A} (right chart in Figure 3), the scatter is dominant, and there is no relationship between *k*_{A} and depth. This is emphasized by the value of the CoD (−0.0671), which is very close to zero. A CoD value equal to zero means that the regression equation does not contribute any further than a model comprised by a single fixed value, with a horizontal line passing through the mean of the *k*_{A} values. This endorses the fact that, with this dataset of 74 wetlands/studies, *k*_{A} was found to be independent from depth. This finding gives support to Kadlec & Wallace (2009) proposal, that *k*_{A} should be the most indicated coefficient for HF wetlands design, because its relationship with depth is null or small.

In the Methods section, it was mentioned that for the second dataset (*n* = 33), based on BOD, there was no information on the length and width of the wetlands, and so the apparent number of tanks in series for the TIS model could not be estimated by Equation (7), and a fixed value of *N* was therefore adopted (*N* = 8). In order to investigate the impact of the utilization of this value of *N* for this second dataset, a sensitivity analysis was made, this time using *N* = 3, for similarity with the *P* = 3 value proposed by Kadlec & Wallace (2009) for the *P–k–C** approach in the estimation of BOD. If *N* had been adopted as 3 in the second dataset (instead of *N* = 8), to match with *P* = 3 adopted in the *P–k–C** approach, the regression equation would still be very similar (*k*_{V} = 0.1701·*h*^{−1.104}), and the exponent ‘*b*’ would be only 4% higher than the one shown in Equation (8). Therefore, the impact of the assumption related to *N* was considered minimal.

Equations (8) and (9) can be indicative of important trends related to the influence of depth on the removal rate coefficients, supported by a considerable number of wetlands/studies. However, the authors feel that they have not reached the status of design equations, especially the one relative to *k*_{V} (Equation (8)), because of the wide scatter in the data and the recognition that more datapoints could be aggregated, in the future, based on additional investigations, in order to make the regression equation even more robust. Ideally, separate regression equations for each type of influent could be developed in the future, as more data are incorporated. Nevertheless, the authors judge that the main trend found here, of *k*_{V} decreasing with depth, is conceptually important and goes to the core of what was intended for the present study.

### Interpretation of the mathematical relationships between depth and removal rate coefficients *k*_{V} and *k*_{A}

*k*

_{V}and

*k*

_{A}as a function of depth

*h*. If one analyzes the structure of the power function of both regression equations, both can be represented by

Based on the volumetric removal rate approach, doubling the wetland depth *h* will double its volume and theoretical hydraulic retention time (*τ*). However, according to Equation (10), if ‘*b*’ is exactly equal to 1.0, then if *h* is doubled, *k*_{V} will be halved [*k*_{V} = *a*·(2 *h*)^{−1} = (*a*/2) (*h*)^{−1}], the product *k*_{V}·*τ* will remain the same [(*k*_{V}/2)·(2*τ*)], and thus, for a steady-state first-order reaction, removal efficiency and effluent concentration will also remain unaffected (see Equation (1)). In regression equation (8), based on the 74 HF wetlands, *b* = 1.059, that is, very close to 1.0. Therefore, in practical terms, one could accept that doubling the wetland depth, while keeping the same surface area, will have a small impact on removal efficiency.

Now, thinking in terms of the areal removal rate approach, if ‘*b*’ is exactly equal to 1.0, then, from Equation (14), *d* = 1 − *b* = 1 − 1 = 0. From Equation (11), one has *k*_{A} = *c* · *h*^{0}, or *k*_{A} = *c*, and finally *k*_{A} = *a* · *ε* (based on Equation (13)). In this theoretical case, depth will have no influence on *k*_{A}, since *k*_{A} will have a fixed value, equal to the regression coefficient ‘*a*’ times the porosity *ε*. If *h* has no influence on *k*_{A}, and since the applied hydraulic loading rate (*q*) is unchanged, because only depth has been changed and the surface area remained the same, then the quotient *k*_{A}/*q* will be the same for any depth, and so will be the expected removal efficiency and effluent concentration (see Equation (2)). From regression equation (8), based on the 74 HF wetlands, *b* = 1.059, and thus *d* = 1.059–1 = 0.059, according to Equation (14). As explained here, this is exactly the exponent ‘*d*’ in regression equation (9) (0.059). This value is very close to zero, indicating that, from the regression analysis based on the dataset of 74 HF wetlands, depth had a very small influence on *k*_{A}. Also in regression equation (9), the coefficient ‘*c*’ is equal to *a* · *ε*, as shown in Equation (13). Therefore, in regression equation (8), the coefficient ‘*c*’ is equal to *a* · *ε* = 0.1676 × 0.35 = 0.0587, which is exactly the ‘*d*’ coefficient in Equation (9).

The practical implications of these intricate but simple mathematical manipulations are made clearer in the next section, which presents a summary graphical representation of the influence of depth and surface area on the removal efficiency of HF wetlands.

### Summary of the influence of depth and surface area on HF treatment wetland performance

Table 3 presents a schematic comparison of a possible search by a designer for strategies for improving performance by increasing the volume of an HF wetland. The table compares a base scenario with three alternatives for doubling the wetland volume: (i) doubling depth and keeping the same surface area; (ii) doubling the surface area by doubling the length, while keeping the same depth and (iii) doubling the surface area by doubling the width, while keeping the same depth. For all scenarios, equations are presented for calculating wetland volume (*V*), surface area (*A*), theoretical hydraulic retention time (*τ*), areal hydraulic loading rate (*q*), volumetric first-order removal rate coefficient (*k*_{V}), areal first-order removal rate coefficient (*k*_{A}), dimensionless performance variables (*k*_{V}·*τ* and *k*_{A}/*q*) and cross-sectional organic loading rate OLR_{CS}.

Variable . | Base scenario . | Comparison with base scenario . | ||
---|---|---|---|---|

Volume of the wetland V (m^{3}) | ||||

Surface area of the wetland A (m^{2}) | ||||

Theoretical hydraulic retention time τ (day) | ||||

Areal hydraulic loading rate q (m/day) | ||||

Volumetric first-order removal rate coefficient k_{V} (day^{−1}) | (decreases with higher h) | (≅ the same, since h is the same) | (≅ the same, since h is the same) | |

Areal first-order removal rate coefficient k_{A} (m/day) | (k_{A} ≅ same; k_{V} decreases; h increases) | (k_{A} ≅ the same; kV ≅ the same and h = the same) | (k_{A} ≅ the same; k_{V} ≅ the same and h = the same) | |

Dimensionless variable k_{V}·τ | (k_{V} decreases; τ increases) | (k_{V} ≅ the same; τ doubles) | (k_{V} ≅ the same; τ doubles) | |

Dimensionless variable k_{A}/q | (q = the same and k_{A} ≅ the same) | (k_{A} ≅ the same; q halves) | (k_{A} ≅ the same; q halves) | |

Cross-sectional organic loading rate OLR_{CS} [g/(m^{2} · day)] | (OLR_{CS} halves because h doubles) | (OLR_{CS} the same) | (OLR_{CS} halves because W doubles) |

Variable . | Base scenario . | Comparison with base scenario . | ||
---|---|---|---|---|

Volume of the wetland V (m^{3}) | ||||

Surface area of the wetland A (m^{2}) | ||||

Theoretical hydraulic retention time τ (day) | ||||

Areal hydraulic loading rate q (m/day) | ||||

Volumetric first-order removal rate coefficient k_{V} (day^{−1}) | (decreases with higher h) | (≅ the same, since h is the same) | (≅ the same, since h is the same) | |

Areal first-order removal rate coefficient k_{A} (m/day) | (k_{A} ≅ same; k_{V} decreases; h increases) | (k_{A} ≅ the same; kV ≅ the same and h = the same) | (k_{A} ≅ the same; k_{V} ≅ the same and h = the same) | |

Dimensionless variable k_{V}·τ | (k_{V} decreases; τ increases) | (k_{V} ≅ the same; τ doubles) | (k_{V} ≅ the same; τ doubles) | |

Dimensionless variable k_{A}/q | (q = the same and k_{A} ≅ the same) | (k_{A} ≅ the same; q halves) | (k_{A} ≅ the same; q halves) | |

Cross-sectional organic loading rate OLR_{CS} [g/(m^{2} · day)] | (OLR_{CS} halves because h doubles) | (OLR_{CS} the same) | (OLR_{CS} halves because W doubles) |

Since the regression equations derived here are principally meant to indicate trends, and not exact relationships, the table presents the influences of the changes with expressions such as ‘increases’, ‘decreases’ or ‘remains approximately the same’, when reporting the impacts on first-order removal rate coefficients (*k*_{V} and *k*_{A}) and removal efficiencies (as given by *k*_{V}·*τ* and *k*_{A}/*q*). For instance, although the regression analyses with the 74 HF wetlands indicated that *k*_{V} would be reduced to approximately half the value if the depth doubled, in the summary table, it is only indicated that *k*_{V} would decrease, without specifying by how much, in order to keep the conclusions more general, and also to take into consideration the uncertainties related to the large scatter in the data and operational conditions of the HF wetlands investigated here.

The table clearly shows that, if one wishes to obtain a higher removal efficiency by increasing the HF wetland volume, this increase should preferably be done by enlarging the surface area, since *k*_{V} and *k*_{A} would remain approximately the same in comparison with the base scenario, and *τ* would increase and *q* would decrease. However, one should not be sure that, in the actual performance, everything will happen as expected *a priori* during the conceptual studies. Hydraulic implications in the geometric shape of a biological reactor and in its inlet and outlet structures are important, and some gains may be offset by occasional hydraulic imperfections that may occur. Clogging may cause disturbances in the flow, possibly inducing dead zones and hydraulic short circuits. Also, microbial communities may be influenced by distinct load distributions in the reactor volume. Therefore, the designer needs to pay attention to additional details and use best practices of design guidelines to conceive the best possible biological reactor. Still, the authors consider that the results of this investigation cast an important light on factors that are decided at the design stage and that have implications on treatment performance.

At the design stage, when looking for alternatives to improve an HF wetland performance, the designer needs to take into account additional practical aspects, exemplified in the three scenarios evaluated in Table 3. In all scenarios (‘*a*’, ‘*b*’ and ‘*c*’), the increase of volume is costly (earth moving and filter media), and this increase is not associated with a proportional improvement in removal efficiency and effluent concentration, as can be seen from Equation (1). Also, in scenarios ‘*b*’ and ‘*c*’, the increase in surface area may be, not only costly but also may bring difficulties for the site allocation of the treatment plant, due to the already high land requirements of an extensive system such as an HF wetland.

### Other practical aspects in the selection of *k*_{V} and *k*_{A}

Another practical question may come to the designer's mind: if the reduction of volume by decreasing depth is not likely to decrease removal efficiency, why not adopt very small HF wetland depths, in order to save in construction costs? This brings another important factor, not dealt with in this study, related to the concept of cross-sectional organic loading rate (refer to the last row of Table 3) and the tendency for clogging. As mentioned in Wallace & Knight (2006) and in Kadlec & Wallace (2009), there is an empirical relationship between the organic loading applied to the inlet cross section of an HF wetland and the potential for problematic accumulation of solids. Initial recommendations were that cross-sectional inlet loadings should be less than 250 g/m^{2}·day BOD for short-term loadings, and Wallace (2014) estimated that less than 100 g/m^{2}·day of BOD would be a safer criterion for long-term loadings. Therefore, if one decreases depth *h*, while keeping the same width *W*, the cross-sectional area (*W*·*h*) will become smaller, and so the cross-sectional organic loading rate will increase, what may bring difficulties for the compliance with the maximum design values in terms of g BOD/m^{2}·day, as discussed earlier.

The balance of positive and negative aspects associated with the possibilities of decreasing or increasing HF wetland depth is complex and, this is one of the reasons that, based on experience, the adoption of liquid depths in the vicinity of 0.50 m are being currently recommended for design purposes (Dotro *et al.* 2017). Although it has been mentioned that Equations (8) and (9) have not reached the status of design criteria, it is instructive to see what the resulting values of *k*_{V} and *k*_{A} for BOD or COD removal could be if one applied both equations for a depth of 0.50 m: *k*_{V} = 0.1676 × 0.50^{−1.059} = 0.35 day^{−1}; *k*_{A} = 0.0587 × 0.50^{−0.059} = 0.061 m/day. As expected, the resulting areal rate coefficient *k*_{A} is within the typical range presented by von Sperling *et al.* (2023) and very similar to the ‘*c*’ coefficient in the regression equation (0.0587), since *k*_{A} was very little influenced by depth.

Another point is worth raising in the comparison between the volumetric and areal removal rate coefficients. The perspective here is for design estimations on effluent concentrations, and in this case the theoretical hydraulic retention time (*τ*) and applied hydraulic loading rate (*q*) are used in Equations (1) and (2). However, if one is studying an existing HF wetland that has been subjected to tracer studies, the mean tracer hydraulic retention time could be used in replacement of the theoretical one in Equations (1) and (6). The mean tracer retention time is a *resultant* measured value that is associated with the effective volume of the reactor, which is usually lower than the total liquid volume due to dead zones and hydraulic short circuits. As a consequence, mean tracer times are usually lower than *τ*. On the other hand, areal-based removal rate coefficients are related to an *applied* areal hydraulic loading rate (*q*). In the representation of this existing wetland that had tracer tests performed, to keep consistency between the volumetric and the areal approaches, in Equation (2) one would need to also correct the areal hydraulic loading rate for the effective area (flow divided by effective area and not flow divided by total available surface area) and assume that the effective area is proportional to the effective volume determined in the tracer tests.

## CONCLUSIONS

The results from the regression analyses of depth versus volumetric (*k*_{V}) and areal (*k*_{A}) first-order removal rate coefficients for the TIS model, using data from HSSF wetlands (*n* = 74) treating domestic wastewater, endorsed the applicability of the areal-based approach for estimating effluent concentrations and removal efficiencies for design purposes. This approach has been largely embraced by the treatment wetland literature, but broad supporting evidence was still scarce before this study.

Based on the results of organic matter (BOD and COD) removal from the overall dataset of 74 wetlands/studies, the areal removal rate coefficient (*k*_{A}) was not heavily influenced by the submerged depth, which is a supporting argument for the adoption of *k*_{A} for the estimation of effluent concentrations. On the other hand, the volumetric removal rate coefficient (*k*_{V}) was strongly affected by depth, in the sense that the increase in depth led to a decrease in *k*_{V}. The volumetric approach can still be used, since it is more widely applied in the general wastewater treatment literature, but its utilization should be in conjunction with the specification of the wetland depth being adopted in the design.

Two regression equations were obtained for the overall dataset of 74 wetlands/studies: *k*_{V} = 0.1676·*h*^{−1.059} and *k*_{A} = 0.0587·*h*^{−0.059}. These regression equations indicated that *k*_{V} was approximately related to the inverse of depth (if depth doubles, *k*_{V} reduces to approximately half), since the exponent of the equation was close to unity. On the other hand, *k*_{A} was almost independent from depth, as given by the exponent close to zero. This reflects the general trends obtained, but of course provides no guarantee that this behavior will be reproduced in individual HF wetlands. The investigations obtained data for BOD and COD, and for diverse types of influents (raw sewage and primary, secondary and tertiary effluents), and all of them had a similar qualitative relationship with depth, that is, *k*_{V} decreased with depth. However, due to the insufficient number of data for each of these applications, no quantitative associations are proposed here.

The results and considerations here are only valid for passive HSSF wetlands and should not be applied to aerated horizontal wetlands and vertical wetlands, whose hydrodynamics are completely different. Also, only organic matter (BOD and COD) has been covered by the experimental data, and no inferences can be made for nitrogen, phosphorus, coliforms and other constituents.

## ACKNOWLEDGEMENTS

The authors express their gratitude to Dr Alisson Carraro Borges (Federal University of Viçosa, Brazil) for making the original raw dataset (COD data) of horizontal flow wetlands available. J.N. acknowledges the 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

*Desempenho de Sistema Alagado Construído para Tratamento de Esgoto Doméstico*

*Utilização de Lírio Amarelo (Hemerocallis flava) em Sistemas Alagados Construídos para tratamento de esgoto doméstico*

*Comparação de modelos de degradação de matéria orgânica em sistemas alagados construídos tratando esgoto sanitário*

*Desempenho de leitos cultivados (‘constructed wetland’) para tratamento de esgoto: contribuições para concepção e operação*

*PhD Thesis*

*k*) for horizontal flow treatment wetlands