Application of divided convective-dispersive transport model to simulate variability of conservative transport processes inside a planted horizontal subsurface flow constructed wetland

This paper offers a novel application of our model worked out in Maple environment to help understand the very complex transport processes in horizontal subsurface flow constructed wetland with coarse gravel (HSFCW-C). We made tracer measurements: Inside a constructed wetland, we had 9 sample points, and samples were taken from each point at two depths. Our model is a divided convective-dispersive transport (D-CDT) model which makes a fitted response curve from the sum of two separate CDT curves showing the contributions of the main and side streams. Analytical solutions of CDT curves are inverse Gaussian distribution functions. This model was fitted onto inner points of the measurements to demonstrate that the model gives better fitting to the inner points than the commonly used convective-dispersive transport model. The importance of this new application of the model is that it can resemble transport processes in these constructed wetlands more precisely than the regularly used convective-dispersive transport (CDT) model. The model allows for calculations of velocity and dispersion coefficients. The results showed that this model gave differences of 4–99% (of velocity) and 2–474% (of dispersion coefficient) compared with the CDT model and values were closer to actual hydraulic behavior. The results also demonstrated the main flow path in the system.


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 to conventional treatment systems (Almuktar et al. 2018).
There are two main types of constructed wetland: 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 an unsaturated, non-permanent flow, and in HFSFCWs there is a saturated, non-permanent flow (Wu et al. 2015;Valipour and Ahn 2016). Our experiments and calculations were performed on HFSFCWs only. We investigated HFSFCWs using coarse gravel as filter medium (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 , 1999Rash and Liehr 1999), evapotranspiration (Galvão et al. 2010;Beebe et al. 2014), and snow melting that 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). Finally, there are the clogging processes 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 andSperling 2013,Tang et al. 2017) .
Due to the factors mentioned above, the hydrodynamic modeling of SFCWs is a challenging task for experts. In these constructions, biofilm activity and root density can be very intense, and more importantly, 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 micro-porous 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 analyze the hydraulic behavior of constructed wetlands (Levenspiel 1972). Scientists have frequently analyzed 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. Tracer tests were taken from each CW. They were filled with different, homogeneously mixed media, gravelly sand and sandy gravel, and both filter materials contained fractions of clay and slit. Samples were taken inside the HSFCW and at the effluent point as well. He concluded that the hydraulic behavior varied considerably within the system. A disadvantageous length to width ratio may have caused this problem as this system was characterized mainly by plug flow with little longitudinal dispersion. The tracer test results showed the main flow paths. The two sides of the CWs were the preferred transport paths in the influent region; in the effluent region, the main flow path was at the middle of the CW. The results presented that both soil filters resulted in heterogeneous flow. Muñoz et al. (2006) performed four tracer studies with bromide as the tracer. They measured four horizontal subsurface flow constructed wetlands. Their results showed that the theoretical retention time was higher in each tracer study than the mean retention time. Tang et al. (2017) studied the hydraulic performance of HSSFCWs. Hydraulic behavior was determined by tracer tests; the tracer was bromide ion. They performed measurements on four CWs. The porosity was between 41 and 44%. The results showed that the theoretical retention times were higher in all tracer tests than the mean retention time. 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 taken at three depths. The results showed that the preferred flow along the bottom layer was with 65-70% of mass flowing along the bottom and 14-18% and 16-17% of mass at the middle and top levels. In this case, the actual residence time was lower than the theoretical time. Richter et al. (2003) measured HSFCW and were the first in the literature who found that the actual HRT could be higher than the theoretical HRT; however, they didn't find the reason for this phenomenon. Bonner et al. (2017) completed a tracer test on laboratoryscale HSFCWs with red fluorescent dye used as tracer. There were 13 sampling points in the CW. The results indicated that the actual residence time was bigger than the theoretical one and that there was some back mixing within the system. The porosity was 43%. We had the same findings based on our transport model (Dittrich and Klincsik 2015a). In this paper we intended to prove that this principle is valid with regard to the inner points of HSSFCWs. Liu et al. (2018) investigated the effect of solids accumulation and root growth on the hydrodynamics of the HSFCWs. They used three laboratory-scale HSFCWs. The tracer was fluorescein sodium. The samples were taken at two points and at three different substrate depths. Their results indicated that the presence of plant root restricted the water flow within the top layer which resulted in the preferential bottom flow phenomenon. The results showed dispersion numbers between 0,09 and 0,16; these numbers were within the acceptable range. Batchelor and Loots (1996) attempted to fit completely stirred series tank reactor (CSTR) and convection-dispersion transport (CDT) models to their tracer test results which yielded bad fitting results, the reason of which the authors could not exactly clarify. Chazarenc et al. (2003) investigated with fitting CSTR and CDT models as well. Their results showed good fittings with CSTR models 9 out of 10 times. Nonetheless, important parameters, for example, porosity and hydraulic conductivity, were estimated values only. King et al. (1997) conducted a conservative tracer analysis of a gravelfilled HSFCW. They fitted CSTR and CDT models as well; bad fittings were found. 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 improvement.
We have developed a new transport model with better fitting properties than the CSTR and CDT models (Dittrich and Klincsik 2015b). This model is a very effective method to obtain more detailed results about the hydraulic behavior of SFCWs. In the present article, our objective is to show that our model can provide an adequate description of complex transport processes inside HSFCW-Cs.

Materials and methods
The tracer measurements were made at a HSFCW-C in Hódmezővásárhely, Hungary. The treatment plant treats 1-1.5 m 3 /day of wastewater from a milk room. The main elements of the technology include a septic tank, the pump system, VHSFW, HSFCW-C, and a polishing pond in the sequence listed.
Scientists have used different tracers: in two studies NaBr (Netter 1994;Tanner and Sukias 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 conservative tracer. The absorption capacity of the filter media for LiCl was tested in the Environmental Technological Laboratory of University of Pécs. The findings indicated that LiCl as a conservative tracer is applicable in the examined construction. For further details about the treatment plant and the tracer tests, consult Dittrich and Klincsik (2015a).
Inside the CW, there were 9 sample points, and the samples were collected at the effluent. These points are demonstrated in Fig. 1. UNICAM Solaar M atomic absorption device was used for the measurement of the LiCl concentration values.
The measured C-t value pairs and other essential measured parameters are summarized in Appendix 1 (Tables 8,9,10,and 11). Four different measurements were performed in time. The measurements obtained S/1, S/2, S/3, and S/4 reference  numbers for easier documentation as shown below. The main data of our own tracer measurements are summarized in Appendix 1 (Table 7).
Our aim was to design a more accurate process than the conventional CDT and CSTR models for simulating transport processes in HSFCW-Cs. We developed an accurate process with the aim to fit the Fréchet distribution function onto effluent tracer test results from HSFCW-C (Dittrich and Klincsik 2015a). Although the Fréchet distribution has good fittings (Dittrich et al. 2020), this method has the disadvantage that the dispersion coefficient cannot be calculated from the fitted Fréchet distribution curve. Therefore, our aim was to develop another method that would provide a solution.
The most popular method for fitting curves to measured response data is the convection-dispersion transport (CDT) model. The 1D equation of the CDT model is where D x [m 2 /h] is the longitudinal dispersion coefficient, q [1/h] is the specific hydraulic loading rate, and ε [-] is the porosity at the time of analysis. The analytical solution of Eq. 1 is with the presumption of a Dirac impulse function at the dose point of the tracer test (Kovács et al. 2002): A three-parameter inverse Gaussian distribution function is shown below: where a and b are the parameters of the distribution function. The mean and standard deviation of the three-parameter inverse Gaussian distribution function are as follows: In this function, parameter c represents delay, so it can substitute for the R parameter that is mathematically separable. With this solution, a parameter is "exempt" from R, so the mean is supplemented: From parameter c, the approximate value of R can be calculated as follows: From Eqs. 6, 7, and 8, the v x , D x , and R parameters can be calculated directly if we can find a good fit for the threeparameter inverse Gaussian distribution curve. This method will be used below for the calculation of transport parameters.
Our model was developed in MAPLE environment; this model fits an inverse Gaussian (IG) distribution function (1st IG curve) onto the main stream and fits the second IG curve onto the side stream while optimizing the main and side stream ratios. For this procedure, we defined parameter s; this value shows the ratio of the first IG curve and the D-CDT curve. Our mathematical method that fits a Fréchet distribution function onto the tracer test curve (Dittrich and Klincsik 2015b) was unified in order to obtain a mathematically stable and fast process.
The fitting criteria for the model were as follows: & At the beginning of the measurement, the tracer concentration is zero all over the CW: & The delay (c 1 ) of the first IG curve is equal to the best fitted Fréchet distribution function. & The delay (c 2 ) of the second IG curve is equal to or higher than the delay of the first IG curve. & The sum of the first and second IG curves has to give the best fitting function (divided convective-dispersive transport (D-CDT) function) (R 2 → 1). & The area below the D-CDT function has to be equal to the area of 100% tracer return.
Our model is described in more detail in the article of Dittrich and Klincsik (2015b).

Results and analysis
In this article, the model presented (Dittrich and Klincsik 2015b) was fitted onto inner points (Fig. 2, points I-IX) of all measurements (S/1, S/2, S/3, S/4). Samples were taken at each inner point at top and bottom levels as we assumed that the transport processes at the top and the bottom would differ. Our aim was to prove that the model gives good fittings at the inner points as well and that our hypothesis is also valid for inner points. Our results provide a deeper understanding of transport processes inside the system and give us a more detailed insight into hydrodynamics of CWs.
All fitted curves had better fitting features than the very well-fitting Fréchet distribution. Even the best curve of the Fréchet distribution had less favorable fitting features than those of our model. Table 1 shows R 2 values for the three model types.
We chose fitting figures that demonstrate considerable differences among the three functions. Figures 2,3,4,5,6,7,8,9,and 10 show the fitting results of the D-CDT model compared with the Fréchet distribution and a conventional CDT model. All images of fitting results are shown in Appendix 3 (Figs. 12,13,14,15,16 and 17).
The red curves on Figs. 2,3,4,5,6,7,8,9,and 10 show the very nice fitting properties of the D-CDT model. The effect of the main stream from the blue curves and the effect of the side streams from black curves can be seen.
Using Eqs. (6, 7 and 8), the v x , D and R of the 1st and 2nd IG curves were calculated. Subsequently, the average R and v x values were calculated and weighted by the areas of IG curves. Calculated results are summarized in Tables 3 and 4. The average velocity means the theoretical value (homogeneous flow regime). The fitting parameters of the measurement S/1-S/4 are summarized in Appendix 2 (Tables 12,13,14,15,16,17,18,and 19).

Discussion
Fitting results Table 1 demonstrates that the D-CDT model gives better fitting than the Fréchet distribution or the CDT model at all points. All average R 2 values were above 0.96. The S/4 top measurement results showed the biggest differences in R 2 values of the D-CDT model (0.97) and the CDT model (0.88).
The S/1 measurement (Table 1) gave the best fitting results of all 4 measurements. This is mainly due to the fact that there was no hydraulic distortion in the CW, as neither the roots nor the formation of the biofilm had any impact upon flow processes. The results show that the D-CDT model has better fittings than the Fréchet distribution and even better than the conventional CDT model fittings.
Upon the S/2 measurement (Table 1), the CW was 1-month old as the results also demonstrate. Even so, almost all R 2 results yielded by the D-CDT model were higher than 0.95.  There was an extreme hydraulic short circuit at the fourth bottom measuring point resulting in a very bad fitting; nevertheless, the D-CDT model provided a better fitting than the CDT model and the Fréchet distribution. As a result of this bad fitting, our model needs to be further developed and extended with a short circuiting part.
Upon the S/3 measurement (Table 1), similar fitting results were obtained (D-CDT model) as upon the second, but the results of the CDT model and Fréchet distribution were considerably worse. For the third, fourth, and sixth bottom measuring points, there were extreme hydraulic short circuits that resulted in very bad fittings. Nonetheless, the D-CDT model Table 3 Calculated transport parameters at S/1 and S/2 top and bottom points based on the fitted D-CDT model proved to provide much better fitting than the CDT model or the Fréchet distribution. Upon the fourth measurement (Table 1), the results of the three functions showed differences. The D-CDT model gave better fittings than the other two methods. This was especially true at the first, third, and seventh top points, where the CDT model showed very bad fittings and the D-CDT model had especially good fittings.

Macro-and micro-porous ratio
When comparing the s values (Table 2), the top value was similar at S/1 and S/2 measurements, while the value of s was reduced at S/3 and S/4 measurements. Overall, it can be said that the rate of micro-porous system volume to total porous volume grew over time from 38 to 57%.  The s values of the bottom measurements presented a gradual decrease until the S/4 measurement. In general, it can also be observed that with regard to the bottom measuring points, the rate of micro-porous system volume to the total porous volume grew over time from 32,9 to 49%.

The role of the main stream and side stream
At the S/1 top and bottom points (Table 3), the side stream had a higher dispersion coefficient than the main stream. It showed significant back mixing in this HSFCW. At points IV, V, VI, VII, and IX, the weighted velocity value was higher than the bottom values which is possible because the CW was 1 day old and the roots and biofilm activity did not affect transport processes.
In Table 3, at the S/2 top and bottom points, the side stream showed a higher dispersion coefficient than the main stream, same as at S/1. At points I, III, IV, V, and VII, the weighted velocity value was higher than the bottom values. At points II, VI, VIII, and IX, these values were the opposite. It means that the roots and the biofilm had been growing and were consequently affecting transport processes. This was possible as the CW was 1 month old. In the cross-section, the middle points (II, V, VIII) showed the lowest velocity values, which led us to the conclusion that it was the root system and biofilm that had the greatest effect on the flow in the middle section from which it follows that the main flow path was on the two sides.
Upon the S/3 measurement showed in Table 4, the CW was 5 months old. The results present that at both the top and bottom points, the side stream had a higher dispersion coefficient than the main stream. At bottom points I, III, VI, VII, and IX, the weighted velocity values were higher than the top values. In the cross-section, the points II, VI, and IX had the lowest velocity values; this means that the root system and biofilm had the most significant impact upon the flow at these points. Consequently, the main flow path was at points III, VI, and VII.
The results gained upon the last measurement (the CW was 9 months old) are shown in Table 4. At both the top and bottom points, the side stream had a higher dispersion coefficient than the main stream, similar to S/1. At points I, II, and III, the top weighted velocity values were higher than the bottom values; therefore, the main flow path was in the first section on the top; at points IV-IX, this was the opposite. This phenomenon is due to the fact that the root system on the top was located much higher and the biofilm activity was much greater than at the bottom. Similar to the results observed by Liu et al. (2018), this resulted in the so-called bottom flow effect, meaning that the main flow path was at the bottom and at the top there was some back mixing, short circuits, and dead zones. Upon all measurements, the preferred transport path is at the influent and effluent regions of the CWs. The reasons for significant back mixing in this CW were as follows: & The distribution pipe was inadequately positioned. The pipe was approximately 0.8 m from the planned place in the direction of the flow, so wastewater could flow back under the distribution pipe. & The seated Carex acutiformis has a globular root system (see Fig. 11); back mixing zones can develop behind these well-defined root zones. This phenomenon is caused by the smaller hydraulic conductivity of fields with high root density, as in such cases the wastewater needs to change flow direction in the filter media.
Our results showed that upon all measurements, the theoretical velocity was higher than the weighted average. This means that the actual HRT was higher than the theoretical HRT (Dittrich and Klincsik 2015b), a result that contradicts results in international literature (Schierup et al. 1990;Breen and Chick's 1995;Tanner and Sukias 1995;King et al. 1997). The values of R are between 1.00-2.67, which show the extent of dead zones in the systems. Table 5 shows the comparison of the main stream D-CDT and the conventional CDT model based on S/1 top and bottom point measurements. At all points, the D-CDT model provided lower porous velocity and higher dispersion coefficient values than the CDT model, with significant differences in some cases. The velocities were 8-71% lower, and dispersion coefficients were 0-352% higher.

Comparison of main stream D-CDT model with conventional CDT model
Upon the S/2 measurement, in most cases (except the VI top point, with a difference of only 2%), the results proved that the D-CDT model provided lower porous velocity and higher dispersion coefficient values than the CDT model as shown in Table 5.The velocities are lower by 0-82%, and dispersion coefficients are higher by 11-394%. Table 6 shows the comparison of main stream D-CDT and conventional CDT models based on S/3 top and bottom point measurements. The D-CDT model yielded lower porous velocity and higher dispersion coefficient values than the CDT model in most cases (except IV bottom point which showed the worst fitting results); consequently, our model needs to be improved and developed further. The velocities were reduced by 7-81%, a range similar to those at S/1 and S/2, and dispersion coefficients were higher by 2-343%.
The results of the final measurement are shown in Table 6. At all top and bottom points, the D-CDT model gave lower porous velocity and higher dispersion coefficient values than the CDT model. The velocities were lower by 11-99%, and dispersion coefficients were higher by 11-390%.

Transport processes controlled by convectivedispersive transport
The diffusion from the perspective of transport processes in the micro-porous system of HSFCW-Cs is not significant. In these systems, there are main stream and side stream regimes with side streams having lower porous velocity and dispersion in the micro-porous system than in the main stream. Nonetheless, it is important to emphasize that these parameters are much higher than diffusion values. In our model, the first inverse Gaussian curve (Figs. 2,3,4,5,6,7,8,9,and 10; blue dotted line) clearly demonstrates the convective-dispersive transport caused by the main stream. The second inverse Gaussian curve (Figs. 2,3,4,5,6,7,8,9,and 10;black line) illustrates that the convective-dispersive transport of water slowed down in microporous systems. The second curve shows the impact of back mixing and dead zones in the side stream combined. These phenomena are significant from the viewpoint of transport processes in such systems. According to our observations, the micro-porous volume grew over time (Table 2). Biofilm activity and root growth have a significant effect upon transport processes, especially with the aging of constructed wetlands. These factors determine flow directions, dead zones, and short circuits in these CWs (e.g., see Table 3, points II, V, and VII). Our model presents these phenomena much more accurately than currently used models (Tables 5 and 6). Our model needs further improvement, especially at points with extreme hydraulic short circuits as our model yielded bad fittings. Nonetheless, our fittings were considerably better than those achieved by the CDT model for example (Table 1; S/3 III, IV, and VI bottom points). The separation of back mixing, dead zones, and side streams are currently being investigated with our model.

Conclusions
A divided convective-dispersive transport (D-CDT) model was designed and constructed with the aim to accurately simulate conservative transport processes in planted, horizontal, subsurface flow constructed wetlands filled with coarse gravel (HSFCW-C). This model makes fitted response curves from the sum of two independent CDT curves, while model optimizing the rate between the main stream and the side streams. In our model, the first inverse Gaussian curve clearly demonstrates the convective-dispersive transport caused by the main stream; the second inverse Gaussian curve illustrates that the convective-dispersive transport of water slowed down in a micro-porous system. The second curve integrates the back mixing and effects of dead zones as well. The diffusion is not significant with regard to transport processes in the microporous system of HSFCW-Cs. The analytical solutions for these two CDT curves are inverse Gaussian distribution functions. Fréchet distribution was used for the rapid optimization of the mathematical procedure (Dittrich and Klincsik 2015b).
The results show that the model is not only adequate and relevant at the effluent points (Dittrich and Klincsik 2015b), but good fittings were achieved at the inner points as well. The D-CDT model proved to give better fittings than the conventional CDT model (Table 1). There are some fittings that yielded considerable differences between the two models. The results indicated that the ratio of micro-porous system volume to the total porous volume grew over time, at the top from 38 to 57% and at the bottom from 32.9 to 49% ( Table 2).
The results also showed that in the 9-month-old CW, the main flow path was at the bottom layers; at the top layer, there were some short circuits and dead zones. Calculated velocity and dispersion coefficients upon using the D-CDT model gave differences of 7-99% (of velocity) and 2-394% (of dispersion coefficient) when compared with the conventional CDT model; results approximated real hydraulic behavior more closely (Tables 5-6).
Overall, it can be said that the designed D-CDT model can be applied to inner points and can thus help discover and understand the hydraulic behavior and transport processes in planted horizontal subsurface flow constructed wetlands.
One of our 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 with their research analyses of transport test results. In our opinion, the outlined statistical method can provide an adequate tool 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 software that could be available for wider application. One of the main directions of our further research is to find other areas where similar research could be performed with success.
Funding Open access funding provided by University of Pécs. The project has been supported by the European Union and co-financed by the European Social Fund under grant agreement no. EFOP-3.6.1.-16-2016-00004. The research was financed by the Higher Education Institutional Excellence Programme of the Ministry for Innovation and Technology in Hungary, 2019, within the framework of the 3rd Thematic Programme of the University of Pécs.
Data availability All data generated or analyzed during this study are included in this published article (and its supplementary information files).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval Not applicable.

Consent to participate Not applicable.
Consent to publish Not applicable.