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 (kV) are traditionally used in conjunction with the theoretical hydraulic retention time. Areal removal rate coefficients (kA) 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 kV and kA. 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, kV decreased with depth of the wetland system. Regression analyses between depth and removal rate coefficients were performed, and the equations indicated that kV was approximately related to the inverse of depth, while kA 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 kV must be provided together with the depth being considered.

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

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 BOD5) 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 (kV) 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 (kA) 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 kV·τ (kV 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 (Qin), the increase in τ is obtained by adopting larger reactor volumes (since τ = (V·ε)/Qin, 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 (kV) 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 kV·τ may not double, because the overall kV 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 (kA, expressed in m/day), coupled with the areal hydraulic loading rate (q, expressed in m/day or m3/(m2·day), corresponding to inflow Qin 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 kA/q, instead of kV·τ.

The volumetric- and areal-based approaches are summarized in Equation (1) (based on kV) and Equation (2) (based on kA), 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):
formula
(1)
where Cin is the influent concentration to wetland (mg/L); Cout 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 (NTIS) (dimensionless); kV is the volumetric removal rate coefficient (day−1); τ is the total theoretical hydraulic retention time in the wetland (τ = V × ε/Q = h×ε/q); Q = flow (m3/day); V is the wetland volume (m3); h is the liquid depth (m); ε is the medium porosity (dimensionless);
formula
(2)
where kA is the areal removal rate coefficient (m/day); q is the applied areal hydraulic loading rate [m3/(m2·day)].
The volume-based and the area-based approaches are interlinked. The connections between the two coefficients may be made by Equations (3)–(5), which are based on the well-known relationship between the theoretical hydraulic retention time (τ), liquid volume (V·ε) and flow (Qin): τ = V·ε/Qin = (A·h)·ε/Qin = h·ε/q:
formula
(3)
formula
(4)
formula
(5)

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, kA 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 kA values from both units were similar and, differently from the findings of García et al. (2004), the shallower unit did not have higher kA 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 kV by half, and so the removal efficiency was kept the same, since the dimensionless product kV·τ ended up being the same. On the other hand, kA values remained the same, and this was one of the reasons for advocating the use of kA instead of kV 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.

In order to obtain a large number of datapoints to perform the regression analyses of depth versus kV and kA, 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.

Table 1

Descriptive statistics of the 74 wetlands/studies that formed the dataset used in this investigation

StatisticsFlow Q (m3/day)Area A (m2)Liquid depth h (m)HRT (day)Areal HLR q (m3/m2·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 
StatisticsFlow Q (m3/day)Area A (m2)Liquid depth h (m)HRT (day)Areal HLR q (m3/m2·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).

The TIS model was used for the estimation of the removal rate coefficients. The estimation of the volumetric removal rate coefficient kV was made by Equation (6), which is a rearrangement of Equation (1), making kV explicit, and by having measured average values of Cin and Cout, assumed values for C*, calculated values of t and estimated values of N:
formula
(6)

The calculation of the areal removal rate coefficient kA 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).

There was no tracer data with information on the NTIS 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.
formula
(7)
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 kV and kA 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.

Descriptive statistics of volumetric (kV) and areal (kA) removal rate coefficients

Before entering the study of the influence of depth on the removal rate coefficients, it is worthwhile to analyze the kV and kA 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 kV and kA 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 kV and kA 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: kV and kA 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 kA 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 kA 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.
Table 2

Descriptive statistics of first-order volumetric (kV) and areal (kA) removal rate coefficients for organic matter (BOD and COD) obtained from the overall dataset with 74 HF wetlands/studies, separated by influent type

StatisticsVolumetric removal rate coefficient kV (day−1)
Areal removal rate coefficient kA (m/day)
RawPrimarySecondaryTertiaryOverallRawPrimarySecondaryTertiaryOverall
Number of data (n10 43 17 74 10 43 17 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 
StatisticsVolumetric removal rate coefficient kV (day−1)
Areal removal rate coefficient kA (m/day)
RawPrimarySecondaryTertiaryOverallRawPrimarySecondaryTertiaryOverall
Number of data (n10 43 17 74 10 43 17 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 
Figure 1

Box plot of volumetric (kV) and areal (kA) first-order organic matter removal rate coefficients obtained in the HF wetlands analyzed in this study, grouping the two datasets with BOD and COD constituents, and specifying by influent type.

Figure 1

Box plot of volumetric (kV) and areal (kA) first-order organic matter removal rate coefficients obtained in the HF wetlands analyzed in this study, grouping the two datasets with BOD and COD constituents, and specifying by influent type.

Close modal

Influence of depth on the volumetric removal rate coefficient (kV), separated by dataset (BOD and COD) and influent type

Figure 2 plots the calculated first-order volumetric removal rate coefficients (kV) 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 kV 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 kV. 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).
Figure 2

First-order removal rate coefficients (kV) separated by datasets (41 datapoints with COD in the left-hand chart and 33 datapoints for BOD in the right-hand chart) and influent types (raw, primary, secondary and tertiary), together with trendlines (power function).

Figure 2

First-order removal rate coefficients (kV) separated by datasets (41 datapoints with COD in the left-hand chart and 33 datapoints for BOD in the right-hand chart) and influent types (raw, primary, secondary and tertiary), together with trendlines (power function).

Close modal

The data from 41 wetlands/studies (left chart) had higher kV 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 kV 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 kV 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 kV 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 kV and kA using the overall dataset

In order to have a broader approach that could sustain a mathematical interpretation of the relationship between saturated depth (h) with kV and kA, 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·xb, 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 kV and kA, as discussed in the next section.
Figure 3

First-order volumetric (kV, left-hand chart) and areal (kA, right-hand chart) removal rate coefficients, using the overall dataset, and putting together BOD and COD and all influent types (raw, primary, secondary and tertiary). Regression curves based on a power function are presented.

Figure 3

First-order volumetric (kV, left-hand chart) and areal (kA, right-hand chart) removal rate coefficients, using the overall dataset, and putting together BOD and COD and all influent types (raw, primary, secondary and tertiary). Regression curves based on a power function are presented.

Close modal
The resulting regression equations are:
formula
(8)
formula
(9)
In the regression analysis between h and kV (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 kV decreasing with the increase in the depth is clearly seen.

On the other hand, in the regression analysis between h and kA (right chart in Figure 3), the scatter is dominant, and there is no relationship between kA 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 kA values. This endorses the fact that, with this dataset of 74 wetlands/studies, kA was found to be independent from depth. This finding gives support to Kadlec & Wallace (2009) proposal, that kA 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 (kV = 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 kV (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 kV 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 kV and kA

Besides analyzing the trends represented by the curves associated with Equations (8) and (9), it is important to discuss about the mathematical interpretation of these two regression equations for kV and kA as a function of depth h. If one analyzes the structure of the power function of both regression equations, both can be represented by
formula
(10)
formula
(11)
Since (see Equation (4)), Equations (4) and (10) can be combined as
formula
(12)
Therefore, the regression coefficients ‘c’ and ‘d’ can be expressed as
formula
(13)
formula
(14)

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, kV will be halved [kV = a·(2 h)−1 = (a/2) (h)−1], the product kV·τ will remain the same [(kV/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 kA = c · h0, or kA = c, and finally kA = a · ε (based on Equation (13)). In this theoretical case, depth will have no influence on kA, since kA will have a fixed value, equal to the regression coefficient ‘a’ times the porosity ε. If h has no influence on kA, 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 kA/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 kA. 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 (kV), areal first-order removal rate coefficient (kA), dimensionless performance variables (kV·τ and kA/q) and cross-sectional organic loading rate OLRCS.

Table 3

Schematic summary of the implications of increasing an HF wetland volume by enlarging depth or surface area, in comparison with a base scenario

VariableBase scenarioComparison with base scenario
     
Volume of the wetland V (m3    
Surface area of the wetland A (m2    
Theoretical hydraulic retention time τ (day)     
Areal hydraulic loading rate q (m/day)     
Volumetric first-order removal rate coefficient kV (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 kA (m/day)   (kA ≅ same; kV decreases; h increases)  (kA ≅ the same; kV ≅ the same and h = the same)  (kA ≅ the same; kV ≅ the same and h = the same) 
Dimensionless variable kV·τ   (kV decreases; τ increases)  (kV ≅ the same; τ doubles)  (kV ≅ the same; τ doubles) 
Dimensionless variable kA/q   (q = the same and kA ≅ the same)  (kA ≅ the same; q halves)  (kA ≅ the same; q halves) 
Cross-sectional organic loading rate OLRCS [g/(m2 · day)]   (OLRCS halves because h doubles)  (OLRCS the same)  (OLRCS halves because W doubles) 
VariableBase scenarioComparison with base scenario
     
Volume of the wetland V (m3    
Surface area of the wetland A (m2    
Theoretical hydraulic retention time τ (day)     
Areal hydraulic loading rate q (m/day)     
Volumetric first-order removal rate coefficient kV (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 kA (m/day)   (kA ≅ same; kV decreases; h increases)  (kA ≅ the same; kV ≅ the same and h = the same)  (kA ≅ the same; kV ≅ the same and h = the same) 
Dimensionless variable kV·τ   (kV decreases; τ increases)  (kV ≅ the same; τ doubles)  (kV ≅ the same; τ doubles) 
Dimensionless variable kA/q   (q = the same and kA ≅ the same)  (kA ≅ the same; q halves)  (kA ≅ the same; q halves) 
Cross-sectional organic loading rate OLRCS [g/(m2 · day)]   (OLRCS halves because h doubles)  (OLRCS the same)  (OLRCS 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 (kV and kA) and removal efficiencies (as given by kV·τ and kA/q). For instance, although the regression analyses with the 74 HF wetlands indicated that kV would be reduced to approximately half the value if the depth doubled, in the summary table, it is only indicated that kV 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 kV and kA 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 kV and kA

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/m2·day BOD for short-term loadings, and Wallace (2014) estimated that less than 100 g/m2·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/m2·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 kV and kA for BOD or COD removal could be if one applied both equations for a depth of 0.50 m: kV = 0.1676 × 0.50−1.059 = 0.35 day−1; kA = 0.0587 × 0.50−0.059 = 0.061 m/day. As expected, the resulting areal rate coefficient kA 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 kA 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.

The results from the regression analyses of depth versus volumetric (kV) and areal (kA) 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 (kA) was not heavily influenced by the submerged depth, which is a supporting argument for the adoption of kA for the estimation of effluent concentrations. On the other hand, the volumetric removal rate coefficient (kV) was strongly affected by depth, in the sense that the increase in depth led to a decrease in kV. 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: kV = 0.1676·h−1.059 and kA = 0.0587·h−0.059. These regression equations indicated that kV was approximately related to the inverse of depth (if depth doubles, kV reduces to approximately half), since the exponent of the equation was close to unity. On the other hand, kA 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, kV 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.

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

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

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)
.
Arceivala
S. J.
1981
Wastewater Treatment and Disposal. Engineering and Ecology in Pollution Control
.
Marcel Dekker
,
New York
, p.
892
.
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
.
Agricultural Engineering Research
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
(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)
.
Revista Brasileira de Engenharia Agrícola e Ambiental-Agriambi
15
(
11
),
1186
1192
(in Portuguese)
.
Chiatti
F. C. F.
&
von Sperling
M.
2012
Influence of retention time, number of ponds and pond depth on nitrogen removal in shallow maturation ponds treating UASB reactor effluent
. In
8th World Congress of IWA – International Water Association
,
September 16–20, 2012
,
Busan, Korea
.
Crites
R. W.
,
Middlebrooks
E. J.
,
Bastian
R. K.
&
Reed
S. C.
2014
Natural Wastewater Treatment Systems
.
CRC Press
,
Boca Raton, FL
, p.
525
.
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
.
Eustáquio
V.
Jr.
,
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 Scientiarum. Technology
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
.
Environmental Technology
.
doi:10.1080/09593330.2020.1769741
.
García
J.
,
Aguirre
P.
,
Mujeriego
R.
,
Huang
Y.
,
Oritz
L.
&
Bayona
J. M.
2004
Initial contaminant removal performance factors in horizontal flow reed beds used for treating urban wastewater
.
Water Research
38
,
1669
1678
.
Kadlec
R. H.
&
Wallace
S. D.
2009
Treatment Wetlands
, 2nd edn.
CRC Press
,
Boca Raton, FL
.
Levenspiel
O.
1999
Chemical Reaction Engineering
, 3rd edn.
John Wiley & Sons, Inc
,
USA
.
Libhaber
M.
&
Orozco-Jaramillio
A.
2012
Sustainable Treatment and Reuse of Municipal Wastewater
.
IWA Publishing
,
London
, p.
557
.
Mara
D. D.
2003
Domestic Wastewater Treatment in Developing Countries
.
Earthscan
,
London
, p.
293
.
Metcalf & Eddy
.
2014
Wastewater Engineering: Treatment and Resource Recovery
, 5th edn.
Metcalf & Eddy/AECOM
. New York, NY, USA, p.
2018
.
Nivala
J.
,
Boog
J.
,
Headley
T.
,
Aubron
T.
,
Wallace
S.
,
Brix
H.
,
Mothes
S.
,
van Afferden
M.
&
Muller
R. A.
2019
Side-by-side comparison of 15 pilot-scale conventional and intensified subsurface flow wetlands for treatment of domestic wastewater
.
Science of the Total Environment
658
,
1500
1513
.
Ren
Y. X.
,
Zhang
H.
,
Wang
C.
,
Yang
Y. Z.
,
Qin
Z.
&
Ma
Y.
2011
Effects of the substrate depth on purification performance of a hybrid constructed wetland treating domestic sewage
.
Journal of Environmental Science and Health, Part A
46
(
7
),
777
782
.
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)
,
November 12–16, 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
.
Ecological Engineering
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, SP
, 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
.
International Journal of Environmental and Analytical Chemistry
91
(
7–8
),
693
707
.
von Sperling
M.
&
Mascarenhas
L. C. A. M.
2005
Performance of very shallow ponds treating effluents from an UASB reactor
.
Water Science and Technology
51
(
12
),
83
90
.
von Sperling
M.
,
Wallace
S. D.
&
Nivala
J.
2022
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
,
November 6–10, 2022
,
Lyon, France
.
von Sperling
M.
,
Wallace
S. D.
&
Nivala
J.
2023
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
.
Science of the Total Environment
859
,
160259
.
doi:10.1016/j.scitotenv.2022.160259
.
Wallace
S. D.
,
2014
Reducing wetland area requirements by using intensification strategies
. In:
Proceedings of the 14th IWA Specialist Group Conference on Wetland Systems for Water Pollution Control
, October 4–8. 2014
(
Zhou
Q.
&
Zhai
J.
, eds).
Tonji University, Chongqing University, and IWA
,
Shanghai
,
China
, pp.
54
64
.
Wallace
S. D.
&
Knight
R. L.
2006
Small-scale Constructed Wetland Treatment Systems: Feasibility, Design Criteria, and O&M Requirements
.
Water Environment Research Foundation (WERF)
,
Alexandria, Virginia
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

Supplementary data