Analysis of conservative tracer measurement results inside a planted horizontal subsurface flow constructed wetland filled with coarse gravel using Frechet distribution

We worked out a method in Maple environment to help understand the difficult transport processes in horizontal subsurface flow constructed wetlands filled with coarse gravel (HSFCW-C). With this process, the measured tracer results of the inner points of a HSFCW-C can be fitted more accurately than with the conventionally used distribution functions (Gaussian, Lognormal, Fick (Inverse Gaussian) and Gamma). This research outcome only applies for planted HSFCW-Cs. The outcome of the analysis shows that conventional solutions completely stirred series tank reactor (CSTR) model and convection-dispersion transport (CDT) model do not describe the internal transport processes with sufficient accuracy. This study may help us develop better process descriptions of very complex transport processes in HSFCW-Cs. Our results also revealed that the tracer response curves of planted HSFCW-C conservative inner points can be fitted well with Frechet distribution only if the response curve has one peak.


Introduction
Constructed wetlands (CWs)-also known as treatment wetlands-are engineered systems for wastewater treatment. Constructed wetlands have a very low or zero energy demand; therefore, operation and maintenance costs are significantly reduced compared with conventional treatment systems (Almuktar et al. 2018).
There are two main types of constructed wetlands: freesurface flow systems (FSF-CW) and subsurface flow systems (SSF-CW). SSF-CWs can be further divided according to the direction of the wastewater flow. Wastewater in SSF-CWs runs either horizontally (in HSSF-CWs) or vertically (in VSSF-CWs) towards the filter media. In VSFCWs, there is unsaturated, non-permanent flow, and in HFSFCWs there is saturated non-permanent flow (Wu et al. 2015;Valipour and Ahn 2016). In our experiments and calculations, only HFSFCWs were considered. We investigated HFSFCWs using coarse gravel as filter media (HFSCW-C). Constructed wetlands can treat a wide variety of polluted water, including municipal, domestic, agricultural or industrial wastewaters (Vymazal 2009).
There are important differences between the ideal and the actual flow. One of the reasons is weather conditions, such as rainfall (Kadlec 1997;Kadlec 1999;Rash and Liehr 1999), evapotranspiration (Galvão et al. 2010;Beebe et al. 2014) and snow melting can have a huge impact on the flow within constructed wetlands. Another important factor is the construction of the CW: the differences in porosity and hydraulic conductivity of filter media in volume and over time (Dittrich and Klincsik 2015a;Licciardello et al. 2019), the active volume of the porous system (Goebes and Younger 2004) and the inlet and outlet positions (Alcocer et al. 2012;Wang et al. 2014;Okhravi et al. 2017). The last is the clogging processes, which are caused by solids accumulation (Carballeira et al. 2016;Lancheros et al. 2017;Liu et al. 2019), biofilm development (Button et al. 2015;Aiello et al. 2016;Vymazal 2018;de Matos et al. 2018), and root density and distribution (de Paoli and von Sperling 2013;Tang et al. 2017).
Due to the factors mentioned above, the hydrodynamic modelling of SFCWs is a challenging task for experts. In these constructions, biofilm activity and root density can be very intensive, and more importantly, the biofilm development and root system growth over time may also be significantly more rapid (Samsó and Garcia 2013;Rajabzadeh et al. 2015). These processes can affect the microporous system, hydraulic conductivity and clogging processes as well (Tanner and Sukias 1995). It is quite challenging and often problematic to estimate these processes or even further, to incorporate these factors into a model.
Conservative tracer tests are commonly used to analyse the hydraulic behaviour of constructed wetlands (Levenspiel 1972). Scientists have frequently analysed SFCWs with conservative tracer tests used as experimental tools to gain more detailed information about the internal hydrodynamics of constructed wetlands (Netter 1994;Suliman et al. 2006; Barbagallo et al. 2011;Wang et al. 2014). Our method was also based on tracer tests. Conservative tracer tests allow for calculations of the hydraulic retention time (HRT) and dispersion coefficient (D) of a hydraulic system. Some scientists have also conducted the same tests in HSFCWs with the same goal. Netter (1994) measured two horizontal subsurface flow constructed wetlands. He conducted tracer tests on each CW. They were filled with different, homogeneously mixed media, gravelly sand and sandy gravel, and both filter materials contained fractions of clay and silt. Samples were taken from inside the CWs and at the effluent point as well. The conclusion was that the hydraulic performance varied considerably inside the system due to the detrimental length to width ratio. Initially, there was plug flow with little longitudinal dispersion in this CW. Breen and Chick (1995) completed a more itemised tracer test as they measured tracer concentration values at the bottom and at the top section of the filter media. Similar hydraulic behaviour was observed as described by Netter (1994); however, the authors attributed it to dead zones and hydraulic shortcuts. Liu et al. (2018) investigated the effect of solids accumulation and root growth on the hydrodynamics of HSFCWs. They used three laboratory-scale HSFCWs. The tracer was fluorescein sodium. Samples were taken from two points and three different substrate depths. The results indicated that the presence of plant root restricted the water flow in the top layer, leading to the preferred, bottom-flow phenomenon. Birkigt et al. (2017) investigated the flow and transport processes on a pilot-scale, horizontal subsurface constructed wetland with tracer tests (bromide, deuterium oxide and uranine). There was one sampling point inside the CW; samples were obtained from three depths. The results showed that the preferred flow distribution consisted of 65-70% of mass flowing along the bottom, and 14-18% and 16-17% of mass at the middle and top levels.
The most commonly used SFCW modelling programs have been HYDRUS2D and FITOVERT (Wang et al. 2011;Kumar and Zhao 2011); nevertheless, these softwares also need further development. Batchelor and Loots (1996) tried to fit completely stirred series tank reactor (CSTR) and convection-dispersion transport (CDT) models too to their tracer test results which yielded bad fitting results; the reason of which the authors did not exactly know. Chazarenc et al. (2003) investigated with fitting CSTR and CDT models as well, which fortunately, resulted in good fitting with CSTR models 9 out of 10 times. Nonetheless, the important parameters, for example, porosity and hydraulic conductivity, were estimated values only. King et al. (1997) conducted a conservative tracer analysis of a gravel-filled HSFCW. They fitted CSTR and CDT models as well; they found bad fittings too. Hydrus 2D uses CSTR and CDT models also at the transport module of the software (Langergraber and Simunek 2011;Langergraber et al. 2009;Toscano et al. 2009); results published nonetheless indicate that the module needs further development.
Several international researchers have shown that CDT and CSTR models do not correlate precisely with tracer test results in HFSCWs (Batchelor and Loots 1996;King et al. 1997;Kumar and Zhao 2011). The CDT model uses Inverse Gaussian distribution, and the CSTR model uses Gamma distribution. Taking into consideration the irregularities in previous studies, we tried to find closer correlations among other distribution function types.

Materials and methods
The tracer measurements were made at a HSFCW-C in Hódmezővásárhely, Hungary. Scientists used different tracers, in two cases NaBr (Netter 1994;Tanner and Suikas 1995), in one of the cases tritium (Netter 1994), in another case a special fluorescent substance (eriochrome acid red) (Breen and Chick 1995) and in four cases LiCl (Schierup et al. 1990;Netter 1994;King et al. 1997;Rash and Liehr 1999).
We chose LiCl as a conservative tracer. The absorption capacity of the filter media for LiCl was tested in the Environmental Technological Laboratory of the University of Pécs. The findings indicated that LiCl is applicable as a conservative tracer in the examined construction. More details of the treatment plant and the tracer test can be found in Dittrich and Klincsik (2015a).
Inside the CW, there were 9 sample points, and samples were collected at the effluent. These points are demonstrated in Fig. 1. The LiCl concentration values of the samplings were measured with a UNICAM Solaar M atomic absorption device.
The results gained at the effluent point have already been published (Dittrich and Klincsik 2015a).
The measured concentration-time value pairs and other relevant measurements are summarised in Appendix 1. We have made four separate measurements at different times and in different seasons. The measurements received S/1, S/2, S/3 and S/4 reference numbers for easier documentation. The main data of our own tracer measurements are summarised in Appendix 1.
We found five applicable distribution function types (Fatigue Life, Lognormal, Frechet, Pearson5 and Inverse Gaussian); for detailed analysis, we used EasyFit program. More information on the selection criteria for the functions can be found in Dittrich and Klincsik (2015a). Subsequently, a more accurate and specific fitting method was established in Maple environment to ensure accurate comparison of results for these functions. This mathematical method is able to fit the functions to the measurement values with specifically defined conditions. Further details are found in Dittrich and Klincsik (2015a). The mathematical procedure was published in Dittrich and Klincsik (2015a). Dittrich and Klincsik (2015a) demonstrated that the Frechet distribution is the best-fitting function to effluent point measurement results. The results show that Frechet had the best average R 2 of the effluent measurement point. Only the Pearson5 R 2 value was sufficiently good, nevertheless, lower than Frechet values The present article aims to investigate which is the bestfitting distribution type in inner points.
In tracer test analysis, scientists do not usually measure porosity; instead, they use the porosity value of newly built filter media before starting the operation or they estimate porosity (Schierup et al. 1990;Tanner and Sukias 1995). In our study (Dittrich and Klincsik 2015a), by measuring the porosity of the analysed HSFCW-C, a very precise analysis was performed. Our results show that the effective porosity of the HSFCW-C decreased by more than 50% in the first 6 months as a result of intense biological activity and root growth. These data were used for the analysis of the transport processes. Detailed information about these results can be found in Dittrich and Klincsik (2015a).

Results and discussion
During the course of our work, the following functions were fitted to the data sets in Appendix 1: Fatigue Life, Frechet, Inverse Gauss, LogNormal, and Pearson5. A customised program in the Maple software was applied for the fittings.   (Table 10). Appendix 3 shows all images of the fittings (Tables 11,12,13,14,15,16,17,and 18).
The S/1 measurement is interesting as the sampling data refer to the CW with only 2 days of age (the installation took place on September 01, 2007). Generally, for points I.-III., we obtained good results. These points were characterised by fast-rising, peaked curves. The second segment (IV.-VI.) had wider, flatter curves due to leakage rate deceleration and mixing processes. Regarding points VII.-IX., the measurement results were no longer included in the run of the function; they also contained smaller and larger jumps; thus, the fitting results significantly deteriorated. Function pictures (Figs. 2, 3, and 4.) support this assessment.
For lower point, VII. was the first where we received bad results as shown in Table 1. There were measurement results, where the functions could not fit well; only the Frechet distribution gave a value above 0.95.
The results of the S/2 measurement reflect the evolution of transport processes of the constructed wetland of 1 month age. The first three measuring points produced similarly favourable fittings. However, in the case of the second segment, only point IV. showed flatter functions. For points V.-VI., we got similarly good fittings than at the first section. This observation is probably a consequence of inhomogeneous flow distribution. The inhomogeneous flow distribution means that in this cross section (IV.-VI.), at point IV., the root growth and the biofilm Fitting the functions was the most difficult at the measurement of S/3. For the first two sections (in Appendix 3: Table 12, S/3 measurement I.-VI. top point figures and in  Table 13 S/3 measurement I. to VI. bottom point figures), the picture of all functions demonstrated that the area under the specified function was too small; only for the third section was it identical to the area drawn by the measurement points. The reason for this observation was that the mechanical pre-treatment of the wastewater did not work well leading to significant clogging in the horizontal flow constructed wetland. Mechanical pretreatment is a septic tank which helps solids settling. The clogged filter media have been replaced with a new filter media of the same type; thus, the pre-treatment problem was solved, so that the subsequent measurement results would no longer be affected by strong clogging processes.
The other reason was that the roots of the tufted sedge had sufficiently developed during the first 5 months in the constructed wetland, resulting in further flow distortions. Due to the development of dead zones, intensive biofilm activity and clogging processes, the role of the secondary stream is significant. The presence of the dead zones is mostly indicated by the poor fitting of the Inverse Gaussian function and the elongated tail length of the curve. This is clearly visible on the following figures (Figs. 8,9,and 10). The red arrows on the following figures (Figs. 8,9,10,11,and 13) show the secondary streams (second peaks). Figures 8, 9, 10, and 11 clearly demonstrate that if the response curve has two peaks, neither function fits well enough. At these measuring points, the use of the divided convective-dispersive model plays an important role (Dittrich and Klincsik 2015b). The results of the last measurement S/4 (May 29, 2008) provide a better picture due to the result yielded by using new filter media of the constructed wetland. Comparing these measurements with  . 11), but as soon as the malfunction stopped, good fittings were achieved (Fig. 12).
Compare images of S/3 and S/4 V. fitted results of top point measurements ( Fig. 11 and 12): The results were completely different; nevertheless, we obtained the expected results. After changing the clogged filter media in the constructed wetland to a new one (same type media), the second peak disappeared (in Appendix 3. . However, it also became apparent that it originally tried to fit a similar shape function into the measurement points. There is a functional problem for this particular measurement that has to be mentioned. There was a two-peak curve which revealed worse fittings; consequently, the use of the divided convective-dispersive model was necessitated. This model could not only fit the first peak but the second as well, so it had much better  (Dittrich and Klincsik 2015b). This point can be seen in Fig. 13.
We investigated the order of functions with the age of the constructed wetland. First, we took the average of each function for each measurement time, as shown in the Table 2. The data in the table refer to the measurements were performed; the duration of each experiment was between 8 and 14 h.
Table 2 clearly demonstrates that the fitting was adequate for each function in the new constructed wetland, but with the ageing of the wetland, the functions became more and more difficult to follow at the measurement points due to the flow distortions caused by root growth or biofilm activity (highlights indicate results that do not reach 0.95). Appendix 2 shows R 2 results for each point. The degree of the fittings is better than those achieved by using conventional models. We determined this value (0.95) as we thought that above this not only the fitting was good enough but also that this value was higher than the ones used in international studies. Figure 14 shows the evolution of R 2 as a function of the age of the CW. The hypothesis that the degree of inaccuracy increases with the age of the CW is apparent when applying Inverse Gauss, LogNormal and Fatigue Life (Fig. 14, yellow, purple and green lines); for the other functions, it is completely different (Fig. 14, blue and red lines).
It can be stated that the results of the fitting deteriorated with time passing through the distortion of the flow, and the Frechet distribution only gave a good fitting when the measurement curve only had one peak. Further research is required if a response curve has two peaks. For this type of modelling, we have been the first to use the Frechet distribution and Pearson5. To date, researchers have only used the Inverse Gaussian distribution, and we got similar fitting results as other international studies.
We investigated the fitting results of the top and bottom measuring points; as we assumed, our results have shown that the values of the top and bottom measurement points may differ according to the position of the unsuitably formed dividing line and the root stratification. First, we measured the length and width of a randomly selected root. The planted Carex Elata has a globular root system (see Fig. 15). Back-mixing zones can form behind these insular root zones causing smaller hydraulic conductivity in such areas; thus, the wastewater needs to change flow direction in the filter media.   Table 3 shows that the bottom points are much more balanced by the fitting of each function. In the case of the top points, however, only Frechet and Pearson5 gave a better fit than in the bottom points; the others were weaker, and when applying the Inverse Gauss, the average R 2 of the top points gave a very bad result. When setting the functions' fitting order, Frechet and Pearson5 again ranked the first two and Inverse Gauss ranked fifth. In the top layer, with slower flow and denser roots, and consequently, more dead zones and more intense biofilm activity, these factors are difficult to adapt to functions. This means that the main flow is at the bottom. Bonner et al. (2016) and Liu et al. (2018) came to a similar conclusion from their results.
In Table 4, the second and third columns contain the fitting results and ranking of the effluent point, as published in Dittrich and Klincsik (2015a), while the fourth and fifth columns contain the results and ranking for inner points. It is striking that the order was comparable with the previous measurement results (Table 4, point X.). The two best-fit functions were Frechet and Pearson5; the worst was the Inverse Gauss. The order of LogNormal and Fatigue Life was interchanged; however, when taking a closer look at the results, it appears that the two values are actually very close. The results met our expectations: the results of the inner points' fittings were very similar to the effluent point fitting results.
The results indicate that regarding internal points, the standard deviation of the R 2 average is higher than at the effluent point, and that the internal points gave worse fitting results.
The average difference between the previous results and the internal points was 0.054. It is important to highlight that none of the functions' average R 2 values reached 0.95; it can be stated, therefore, that none of them fitted perfectly with the measuring points.

Conclusions
The purpose of our research was to find better-fitted distribution functions than those conventionally used to our conservative tracer test results at the inside points of a Hungarian HSFCW-C. We fitted 5 distribution functions in the Maple software onto tracer test results of our inner points. These 5 function types were chosen from among a large amount of distribution functions (Dittrich and Klincsik 2015a).
We have determined that the Inverse Gauss function ranked 5th in the order of alignment of the functions. In two cases, it was necessary to modify the parameters manually to fit the specified points. The analytical solution of the CDT model is an Inverse Gaussian distribution function. Therefore, it seems clear that the normal CDT model cannot precisely generate a correctly fitting correlation, as the R 2 values did not reach 0.95, and below this value, the fitting did not give the expected results. The error of the CDT model increases with the age of the CW.
The Fatigue Life and LogNormal distributions in the order of alignment will always be 3-4, which means that these two distributions take the third and fourth places at both the effluent and inner points (Table 3 and Table 4). They can be ignored in further investigations, due to bad average R 2 values.   The first two places were achieved by Frechet and Pearson5. The averages of R 2 and the fitting images of the functions indicate that the Frechet distribution incorporates the measurement points more eloquently than the Pearson5 distribution.
With this process, we proved that the planted HSFCW-C conservative tracer response curves at inner points of CW demonstrate a Frechet distribution. This result is identical to results published about the effluent point of the same CW (Dittrich and Klincsik 2015a). The Frechet distribution proved to be the best fitting only where the measured curve had one peak. Where the measured curve had two peaks, the Frechet distribution did not fit sufficiently well; thus, further research is needed. Investigating the top and bottom measurement points, we found that the fitting results at the top measuring points revealed much worse fittings than the bottom measuring points. This is possible as the top layer is characterised by slower flows, denser roots, more dead zones and more intense biofilm activity.
We carried out similar measurements in another constructed wetland in Pécs, and aim to publish the results in another article.
One of our main goals with this fitting procedure in Maple environment was to provide a novel, adaptable method of analysis for other types of hydraulic regimes and thereby, to aid scientists in their analysis of transport test results. In our opinion, the presented statistical method can be used for a deeper understanding of several hydrodynamic problems for the solution of which traditional methods have not been successful, mainly hydraulic leakage problems in other media. Our further research direction is to develop a general software that would allow a wider application. One of the main directions of our future research is to find other areas where similar research success could be achieved.        (   Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.