Modeling of Pollutants Removal in Subsurface Vertical Flow and Horizontal Flow Constructed Wetlands

Reject water is a by-product of every municipal and agro-industrial wastewater treatment plant (WWTP) applying sewage sludge stabilization. It is usually returned without pre-treatment to the biological part of WWTP, having a negative impact on the nitrogen removal process. The current models of pollutants removal in constructed wetlands concern municipal and industrial wastewater, whereas there is no such model for reject water. In the presented study, the results of treatment of reject water from dairy WWTP in subsurface vertical flow (SS VF) and subsurface horizontal flow (SS HF) beds were presented. During a one-year research period, SS VF bed reached 50.7% efficiency of TN removal and 73.8% of NH4-N, while SS HF bed effectiveness was at 41.4% and 62.0%, respectively. In the case of BOD5 (biochemical oxygen demand), COD (chemical oxygen demand), NH4-N, and TN (total nitrogen), the P-k-C* model was applied. Multi-model nonlinear segmented regression analysis was performed. Final mathematical models with estimates of parameters determining the treatment effectiveness were obtained. Treatment efficiency increased up to the specific temperature, then it was constant. The results obtained in this work suggest that it may be possible to describe pollutant removal behavior using simplified models. In the case of TP (total phosphorus) removal, distribution tests along with a t-test were performed. All models predict better treatment efficiency in SS VF bed, except for TP.


Introduction
In Europe, the biggest municipal and agro-industrial wastewater treatment plants (WWTPs) utilize anaerobic sewage sludge digestion with biogas production [1].The possibility of heat and electric energy production from biogas is an advantage of such sewage sludge treatment, which decreases the cost of waste treatment significantly.
One of the major problems connected with anaerobic digestion chamber exploitation is reject water, which is generated in the process of dewatering stabilized sewage sludge.It is usually returned without pre-treatment to the biological part of a WWTP, having a negative impact on the treatment process and causing problems with reaching high effectiveness of nitrogen removal [2][3][4].Reject water from anaerobic sewage sludge digestion is characterized by a high concentration of nitrogen in the form of ammonium nitrogen (NH 4 + -N) and irregular flow because of periodically working devices for sewage sludge dewatering [5][6][7][8].The problem of reject water treatment in a conventional sludge activated system is caused by low biodegrability, indicated by proportions of easily biodegradable organic matter expressed by biochemical oxygen demand (BOD 5 ) to total organic matter expressed by chemical oxygen demand (COD).Constructed wetlands are considered environmentally friendly technology [9][10][11].They have been used worldwide for the treatment of municipal and industrial wastewater, as well as reject water from anaerobic sewage sludge digestion in municipal WWTPs and landfill leachate, which is also characterized by high NH 4 + -N concentration [12][13][14][15].Reject water generated during anaerobic digestion in dairy WWTPs differs from that from municipal WWTPs.There is a lack of experiments concerning constructed wetlands treatment of reject water generated during anaerobic digestion in dairy WWTPs [16].
The aim of the research was to define parameters, with their expected values, of the mathematical models describing the pollutants removal in subsurface vertical flow (SS-VF) and subsurface horizontal flow (SS-HF) constructed wetlands treating reject water from anaerobic sludge digestion in dairy WWTPs.The P-k-C* model was applied [12,17,18].The issue of modeling the functioning of constructed wetlands is one of the most important challenges in that field [19,20].The current models of pollutant removal concern the municipal and industrial wastewater treatment, whereas there is no such model for reject water treatment, so the necessity for such research goals was confirmed.

Study Sites
The research was conducted over a twelve-month period using a pilot-scale plant built in the biggest dairy plant in Poland, which Person Equivalent (PE) is up to 500,000, located in Wysokie Mazowieckie in the Podlaskie province.The average air temperature is 6.8 • C and precipitation is 562 mm year −1 .
A problem with high concentration of nutrients in the reject water was observed after changing aerobic sewage sludge stabilization to anaerobic.After sewage sludge dewatering with centrifuge, reject water was returned to the main sewage line without separate treatment.In comparison with reject water obtained during aerobic sewage sludge stabilization, the concentration of ammonium nitrogen increased significantly [2].
The research installation (N 52.54.20-E 22.299.20)was designed by the authors using earlier experience with reject water treatment.Its main elements are SS VF and SS HF beds with 5 m 2 surface area and 0.8 m height (Figure 1).The beds were used in parallel.In addition, the installation includes a sedimentation and retention tank, outflow, and sampling points (I, II, III). Figure 2 presents the cross section of SS VF and SS HF beds.working devices for sewage sludge dewatering [5][6][7][8].The problem of reject water treatment in a conventional sludge activated system is caused by low biodegrability, indicated by proportions of easily biodegradable organic matter expressed by biochemical oxygen demand (BOD5) to total organic matter expressed by chemical oxygen demand (COD).Constructed wetlands are considered environmentally friendly technology [9][10][11].They have been used worldwide for the treatment of municipal and industrial wastewater, as well as reject water from anaerobic sewage sludge digestion in municipal WWTPs and landfill leachate, which is also characterized by high NH4 + -N concentration [12][13][14][15].Reject water generated during anaerobic digestion in dairy WWTPs differs from that from municipal WWTPs.There is a lack of experiments concerning constructed wetlands treatment of reject water generated during anaerobic digestion in dairy WWTPs [16].
The aim of the research was to define parameters, with their expected values, of the mathematical models describing the pollutants removal in subsurface vertical flow (SS-VF) and subsurface horizontal flow (SS-HF) constructed wetlands treating reject water from anaerobic sludge digestion in dairy WWTPs.The P-k-C* model was applied [12,17,18].The issue of modeling the functioning of constructed wetlands is one of the most important challenges in that field [19,20].The current models of pollutant removal concern the municipal and industrial wastewater treatment, whereas there is no such model for reject water treatment, so the necessity for such research goals was confirmed.

Study Sites
The research was conducted over a twelve-month period using a pilot-scale plant built in the biggest dairy plant in Poland, which Person Equivalent (PE) is up to 500,000, located in Wysokie Mazowieckie in the Podlaskie province.The average air temperature is 6.8 °C and precipitation is 562 mm year −1 .
A problem with high concentration of nutrients in the reject water was observed after changing aerobic sewage sludge stabilization to anaerobic.After sewage sludge dewatering with centrifuge, reject water was returned to the main sewage line without separate treatment.In comparison with reject water obtained during aerobic sewage sludge stabilization, the concentration of ammonium nitrogen increased significantly [2].
The research installation (N 52.54.20-E 22.299.20)was designed by the authors using earlier experience with reject water treatment.Its main elements are SS VF and SS HF beds with 5 m 2 surface area and 0.8 m height (Figure 1).The beds were used in parallel.In addition, the installation includes a sedimentation and retention tank, outflow, and sampling points (I, II, III). Figure 2 presents the cross section of SS VF and SS HF beds.The SS VF bed was built of three layers of stones 0.2 m (fraction 16-60 mm), gravel 0.4 m (fraction 2-16 mm), and sand 0.2 m (fraction 0.5-2 mm).Two pipes (diameter 50 mm) with perforation 2 mm were used for passive aeration.The SS-HF bed was built of gravel (2-16 mm) and stones (16-60 mm).The gravel hydraulic conductivity was 4×10 −2 ms −1 .Both beds were planted with reeds (Phragmites australis).Reject water was taken directly from centrifuge which dewatered sludge after anaerobic digestion in a dairy WWTP.

Sampling and Analytical Procedures
The study was carried out between April 2015 and March 2016.Both beds were supplied from sedimentation and retention tanks with the same amount of reject water, which allowed comparison of these beds' effectiveness.The daily hydraulic load was 0.1 m 3 m −2 d −1 in both SS VF and SS HF beds.Hydraulic retention time (HRL) for SS HF bed was approximately 8 days.Samples were collected three times a month (influent to SS-VF and SS-HF and effluents from both beds).The air temperature during the research period varied from −11 to 26 °C, while reject water temperature varied from 4 °C to 20 °C.
The basic physical and chemical analyses were performed: BOD5, COD, total organic carbon (TOC), total suspended solids (TSS), total Kjeldahl nitrogen (TKN), ammonium nitrogen (NH4 + -N), nitrate nitrogen (V) (NO3-N), nitrite nitrogen (III) (NO2-N), total phosphorus (TP), dissolved oxygen, and alkalinity.TN value was calculated as a sum of TKN, NO3 -N, and NO2 − -N.To evaluate biodegrability of reject water, BOD5/COD and BOD5/TN ratios were determined.Determinations were conducted in a certified laboratory in accordance with the procedures set out in the Regulation of the Environmental Protection Minister [21] from November 18, 2014, and in accordance with the American Public Health Association (2005) [22].

Modeling of Pollutants Removal
Statistical analysis of NH4 + -N and TN, as well as BOD5; COD removal was performed, based on the P-k-C * model [12]: The SS VF bed was built of three layers of stones 0.2 m (fraction 16-60 mm), gravel 0.4 m (fraction 2-16 mm), and sand 0.2 m (fraction 0.5-2 mm).Two pipes (diameter 50 mm) with perforation 2 mm were used for passive aeration.The SS-HF bed was built of gravel (2-16 mm) and stones (16-60 mm).The gravel hydraulic conductivity was 4 × 10 −2 ms −1 .Both beds were planted with reeds (Phragmites australis).Reject water was taken directly from centrifuge which dewatered sludge after anaerobic digestion in a dairy WWTP.

Sampling and Analytical Procedures
The study was carried out between April 2015 and March 2016.Both beds were supplied from sedimentation and retention tanks with the same amount of reject water, which allowed comparison of these beds' effectiveness.The daily hydraulic load was 0.1 m 3 m −2 d −1 in both SS VF and SS HF beds.Hydraulic retention time (HRL) for SS HF bed was approximately 8 days.Samples were collected three times a month (influent to SS-VF and SS-HF and effluents from both beds).The air temperature during the research period varied from −11 to 26 • C, while reject water temperature varied from 4 • C to 20 • C.
The basic physical and chemical analyses were performed: BOD 5 , COD, total organic carbon (TOC), total suspended solids (TSS), total Kjeldahl nitrogen (TKN), ammonium nitrogen (NH 4 + -N), nitrate nitrogen (V) (NO 3 -N), nitrite nitrogen (III) (NO 2 -N), total phosphorus (TP), dissolved oxygen, and alkalinity.TN value was calculated as a sum of TKN, NO 3 -N, and NO 2 − -N.To evaluate biodegrability of reject water, BOD 5 /COD and BOD 5 /TN ratios were determined.Determinations were conducted in a certified laboratory in accordance with the procedures set out in the Regulation of the Environmental Protection Minister [21] from November 18, 2014, and in accordance with the American Public Health Association (2005) [22].

Modeling of Pollutants Removal
Statistical analysis of NH 4 + -N and TN, as well as BOD 5 ; COD removal was performed, based on the P-k-C * model [12]: Water 2019, 11, 180 where C out -output concentration (g/m 3 ), C in -input concentration (g/m 3 ), C*-background concentration (g/m 3 ), k(T)-chemical reaction coefficient (temperature dependent), q-hydraulic load (m/d), P-number of tanks in series.Model (1) can be further simplified if some assumptions are made.The apparent removal efficiency η app can be defined as: This apparent efficiency model's real efficiency η in limit C * → 0 .If limit P → ∞ is taken, the right side of Equation (1) simplifies: Chemical reaction coefficient k is dependent on reject water temperature, using direct or modified first order Arrhenius dependency: where k d , k m -direct and modified chemical reaction coefficients; K 20 -direct or modified chemical reaction coefficients normalized for 20 • C; θ-temperature coefficient; θ m -modifier for temperature coefficient for temperatures less than specific temperature T k , where trend changes; min(•, •)-minimum function of 2 arguments.Modified first order dependency allows for an additional change of the reaction coefficient below an estimated specific temperature.It was introduced to better fit modeled data, while not changing the basic rule of first order linear dependency on a logarithmic scale.
By using all available combinations of previously described equations, 16 models with different parameters subsets were estimated.The final parameter subset selection was performed using methodology given by Burnham & Anderson [23].A detailed description of fitted models, their selection and fitting procedure are presented in Appendix A.
It is worth noting that a full analysis of a P-k-C* model is difficult and can lead to unexpected results [24].In most of the cases, plug flow models can be too simple to describe the dynamics of particular pollutant removal.Sensitivity analysis of fitted models in literature [25] may suggest that some parameters, like P, may be insensitive and set to values obtained from earlier studies.There are some earlier results suggesting that temperature dependency can change [16].

Treatment Efficiency and Load Removal
Table 1 presents the characteristics of reject water before and after treatment with SS VF and SS HF beds during research period, while Table 2 shows removed pollutants load.
The BOD 5 /COD and BOD 5 /TN ratios give information about biodegradability [26].Its value can be used to assess the reject water susceptibility for high efficiency of biological treatment in conventional systems (e.g., activated sludge method).An average BOD 5 /COD ratio was 0.59, while BOD 5 /TN ratio 0.43.Low BOD 5 /COD ratio pointed out low degradability of the organic compounds.Reject water from dairy WWTP with anaerobic sewage sludge digestion was discovered to have a higher BOD 5 /COD ratio than municipal WWTPs.In municipal, it ranges from 0.25 to 0.32 for reject water from the WWTP in Gdansk [27], and 0.2 for reject water from the WWTP in Minworth [8].The WWTP in Gdansk records BOD 5 /TN ratios at the level of 0.37 up to 0.54.In this case, the conventional path of nitrogen removal cannot be applied as the content of easily biodegradable organics is insufficient [28,29].A very high concentration of ammonium nitrogen should be viewed as the main reason for this.Through the analysis of the effect of organic matter removal (g m −2 d −1 ) expressed by BOD 5 , COD, and TOC (Table 2), similar values to household sewage treatment with constructed wetlands were achieved [10].In the case of TN and NH 4 + -N, higher efficiency was observed mainly due to high concentration of nitrogen in reject water before treatment in SS VF and SS HF beds.The TP removal effect was similar to the observed in the case of household and municipal sewage treatment.

Modeling of Pollutants Removal
Statistical analysis of NH 4 + -N and TN, as well as BOD 5 , COD removal, presented in Appendix B, revealed existence of 3 group models, each with stable core parameters, which dominate the variability of logarithmic likelihood and residuals distribution.The most-preferred parameter set was T k , θ m , K 20 -a model without temperature dependency after T k , with Akaike weight greater than 0.94 in all but other cases except in one (TN removal in SS HF bed).Residual analysis suggests also T k , θ m , K 20 model selection in this case.Selected models are presented graphically in Figure 3, using efficiency scale.
variability of logarithmic likelihood and residuals distribution.The most-preferred parameter set was , , -a model without temperature dependency after , with Akaike weight greater than 0.94 in all but other cases except in one (TN removal in SS HF bed).Residual analysis suggests also , , model selection in this case.Selected models are presented graphically in Figure 3, using efficiency scale.In the cases of BOD5, COD, TN, and NH4 + -N, the treatment efficiency increased until a specific temperature was reached, and then it was constant.Results obtained in this work suggest that it may be possible to effectively describe pollutant removal behavior in a consistent way using simplified models.Any additional parameter, with assumption of proper overall fit, will have a large confidence interval, suggesting overall insensitivity.Simplicity of obtained models does not interfere with their statistical validity.
Created models show that simple ratios obtained over the whole period of the research study are insufficient to properly describe behavior of beds, and a more dedicated approach, like P-k-C* modeling, is necessary.Only ratios in stable periods of bed work (after specific temperature ) are sufficient to describe it.
Values of specific temperature were marked in place of trend change and on the upper and lower section of figures along with 95% confidence intervals (CI).Parameter estimates of the fitted models, along with their 95% confidence intervals, and residual sum of squares were presented in Table 3.In the cases of BOD 5 , COD, TN, and NH 4 + -N, the treatment efficiency increased until a specific temperature T k was reached, and then it was constant.Results obtained in this work suggest that it may be possible to effectively describe pollutant removal behavior in a consistent way using simplified models.Any additional parameter, with assumption of proper overall fit, will have a large confidence interval, suggesting overall insensitivity.Simplicity of obtained models does not interfere with their statistical validity.Created models show that simple ratios obtained over the whole period of the research study are insufficient to properly describe behavior of beds, and a more dedicated approach, like P-k-C* modeling, is necessary.Only ratios in stable periods of bed work (after specific temperature T k ) are sufficient to describe it.
Values of specific temperature T k were marked in place of trend change and on the upper and lower section of figures along with 95% confidence intervals (CI).Parameter estimates of the fitted models, along with their 95% confidence intervals, and residual sum of squares were presented in Table 3.
Figure 4 presents standard a box-whisker plot of TP removal efficiency and reject water temperature during the whole research period.Mean values are described as red dots inside each box, and horizontal line represents median.Dashed lines inside efficiency vs. temperature plot represent mean values.
The calculated treatment efficiency of TP removal was on average 25.42% (SS VF) and 37.20% (SS HF).The efficiency of TP removal in SS VF bed was stable and higher than in SS HF bed.Phosphorous removal occurs mainly in the processes of chemical precipitation and sorption, which do not depend on temperature [12].A distribution analysis of TP treatment efficiency revealed no significant deviation from normality (Shapiro-Wilk normality test, for SS HF: test statistics 0.952, p-value = 0.1, for SS VF: test statistics 0.96, p-value = 0.19).An F-test revealed no statistically significant difference between the variances (F test statistics 1.74, numerator and denominator degrees of freedom-37, p-value = 0.095).T-test revealed statistically significant differences between mean treatment efficiency of TP (t statistics 9.52, 74 degrees of freedom, p-value < 2×10 −14 ).For all performed tests, significance level was set to α = 0.05.The calculated treatment efficiency of TP removal was on average 25.42% (SS VF) and 37.20% (SS HF).The efficiency of TP removal in SS VF bed was stable and higher than in SS HF bed.Phosphorous removal occurs mainly in the processes of chemical precipitation and sorption, which do not depend on temperature [12].A distribution analysis of TP treatment efficiency revealed no significant deviation from normality (Shapiro-Wilk normality test, for SS HF: test statistics 0.952, pvalue = 0.1, for SS VF: test statistics 0.96, p-value = 0. 19).An F-test revealed no statistically significant difference between the variances (F test statistics 1.74, numerator and denominator degrees of freedom-37, p-value = 0.095).T-test revealed statistically significant differences between mean treatment efficiency of TP (t statistics 9.52, 74 degrees of freedom, p-value < 2×10 −14 ).For all performed The mathematical models found practical implementation.An application, which allows determining the effectiveness of reject water treatment in SS VF and SS HF constructed wetlands depending on the reject water temperature, was developed (Figure 5).freedom-37, p-value = 0.095).T-test revealed statistically significant differences between mean treatment efficiency of TP (t statistics 9.52, 74 degrees of freedom, p-value < 2×10 −14 ).For all performed tests, significance level was set to α = 0.05.
The mathematical models found practical implementation.An application, which allows determining the effectiveness of reject water treatment in SS VF and SS HF constructed wetlands depending on the reject water temperature, was developed (Figure 5).

Conclusions
The study showed a high efficiency of SS VF and SS HF beds for removing main pollutants from the reject water generated during anaerobic sewage sludge digestion in dairy WWTP.Results for SS

Conclusions
The study showed a high efficiency of SS VF and SS HF beds for removing main pollutants from the reject water generated during anaerobic sewage sludge digestion in dairy WWTP.Results for SS VF and SS HF, respectively, were: COD 79.79% and 75.34%, for BOD 5 82.12% and 74.13%, for TN 51.46% and 42.58%, for NH 4 + -N 74.06% and 62.52%, for TP 25.42% and 37.20%.A higher efficiency of main organic pollutants removal was observed in the case SS VF bed during the whole research.The efficiency of TP removal was stable and higher for the SS HF bed.
The results have been used for modeling the work of SS VF and SS HF beds for reject water treatment.Models for determining the effectiveness of treatment depending on the temperature and type of beds were presented.P-k-C* model was sufficient to describe the observed variability of constructed wetland behavior.In the cases of BOD 5 , COD, TN, and NH 4 + -N, the treatment efficiency increased to the specific temperature, then it was constant.Results obtained in this work suggest that it may be possible to effectively describe pollutant removal behavior using simplified models.All models predict better treatment efficiency in SS VF bed, except for TP.In the case of TP removal, distribution tests along with a t-test were performed.The created models show that simple ratios obtained over the whole period of study are insufficient to properly describe behavior of beds.A more dedicated approach, like P-k-C* modeling, is necessary.Only ratios in stable periods of bed work (after specific temperature) are sufficient to describe it.
An application in the form of software for modeling the efficiency of reject water treatment from sewage sludge digestion in dairy WWTPs was prepared on the basis of the results from the research.

Figure 1 .
Figure 1.Scheme of research installation, with sampling points (I, II, III).Figure 1. Scheme of research installation, with sampling points (I, II, III).

Figure 1 .
Figure 1.Scheme of research installation, with sampling points (I, II, III).Figure 1. Scheme of research installation, with sampling points (I, II, III).

Figure 2 .
Figure 2. Cross section of SS VF (subsurface vertical flow) and HF (subsurface horizontal flow beds) beds.

Figure 2 .
Figure 2. Cross section of SS VF (subsurface vertical flow) and HF (subsurface horizontal flow beds) beds.

Figure 4 .
Figure 4. TP removal efficiency and reject water temperature during the research period.

Figure 4 .
Figure 4. TP removal efficiency and reject water temperature during the research period.

Figure 5 .
Figure 5.An application "Modeling of reject water treatment efficiency".

Figure 5 .
Figure 5.An application "Modeling of reject water treatment efficiency".

Table 1 .
The characteristics of reject water before and after treatment in SS VF and SS HF beds (38 research series, 114 samples).

Table 2 .
Removed pollutant loads (mean values).After the treatment, the average ratios of BOD 5 /COD were 0.53 for SS VF bed and 0.63 for SS HF bed, and BOD 5 /TN were 0.16 and 0.20, respectively.

Table 3 .
Parameters of selected models for BOD 5 , COD, NH 4 + -N and TN.
Water 2019, 11, x FOR PEER REVIEW 8 of 19

Table A1 .
Fitted models for BOD 5 with their statistics.

Table A2 .
Fitted models for COD with their statistics.

Table A3 .
Fitted models for NH 4 + -N with their statistics.

Table A4 .
Fitted models for TN with their statistics.