Evaluation of Different Methods to Assess the Hydraulic Behavior in Horizontal Treatment Wetlands

While there have been numerous studies on the rate and development of clogging in horizontal subsurface treatment wetlands (HSTWs) and, consequently, the effects on its hydraulic characteristics, research has not shown a clear understanding of the processes. The existing methods for measuring the impact of clogging provide limited information on the extension and degree of the phenomenon. This study aimed to evaluate the capacity of various measurement techniques to assess the degree and variation in space and time of clogging in HSTWs. Hydraulic conductivity at saturation (Ks) measurements were conducted using a newly implemented scheme, the drainage equation method, and traditional tracer tests, which were carried out in a full-scale HSTW system, located in Sicily, Italy, during 2019. After five years of operation, the results highlighted a severe decrease in Ks (<1000 m day−1) in the inlet zone (despite the fact that the filter gravel was replaced in 2017), a very high reduction of Ks along the central path inside the bed, a nonuniform flow through the HSTW, the presence of stagnant zones, and a reduction of the porosity of the bed gravel. Nonetheless, the mean values of the physical–chemical and bacteriological parameters at the hybrid treatment wetland (hybrid TW) outlet indicated that the partial clogging had no significant effect on the quality of the discharged water. Moreover, the results obtained using the different measurement techniques (in terms of both the Ks values and the flow distribution inside the bed) were consistent with each other and with results obtained previously for the same system. Finally, the most efficient combination of methods to assess clogging in HSTWs was identified.


Introduction
Bed clogging of horizontal subsurface treatment wetlands (HSTWs) is inevitable and not a new phenomenon. Its development and impacts on the design and operation of these systems must be taken into account to avoid reductions of their lifetime [1][2][3]. Clogging is also a complex phenomenon, because sets of mechanisms act differently during the lifetime of HSTWs [4]. Cleaning the bed media in the inlet zone of HSTWs (where microbial biomass formation primarily occurs) may be unavoidable and should be regularly scheduled. On the other hand, designers have to take into account that the clean-bed hydraulic conductivity at saturation (K s ) will change significantly once the system is put into operation and that HSTWs must nonetheless be able to operate efficiently. The K s value, which is largest at start-up, will decrease with time as plant roots, microbial biofilms, and chemical precipitates This paper aims (1) to assess the degree and variation in space and time (also by comparison with previous experimental campaigns) of clogging in the HSTW system functioning as a secondary wastewater treatment system of the IKEA ® store, and (2) to identify the most efficient way to evaluate the phenomenon in HSTW systems. In order to reach these objectives, an experimental campaign integrating traditional and innovative methods to assess the hydraulic behavior in HSTWs was carried out in the horizontal bed of the hybrid TW system of the IKEA ® store, which had already been partially clogged in the past. In particular, traditional tracer tests were compared with K s measurements in a newly implemented scheme (i.e., by using a pervious permeameter with a calibrated equation) and the drainage equation method. To ensure that the HSTW conditions were the same, the three methods were carried out consecutively in February 2019; it was not possible to perform them simultaneously due to the characteristics of the individual methods.

Full-Scale HSTWs Characterization
The selected HSTW was chosen as a case study to assess clogging using various measurement techniques. The HSTW system is part of the hybrid TW system that functions as a secondary wastewater treatment system for the IKEA ® store, located in the industrial district of Catania, Italy (37 • 26 54.2" N 15 • 02 05.2" E, 11 m a.s.l., Figure 1). The area is semiarid and characterized by average annual precipitation of about 760 mm. The values for temperature and humidity were in the range of 1.6-40 • C and 8-100%, respectively, in 2018. The average monthly precipitation was 55 mm during January and February 2019. No rains occurred in the days during which the measurements were made or on the three previous days. The treatment plant includes a screening unit and a sequential batch reactor (SBR) as the primary treatment. The SBR was designed for treating wastewater produced by toilets and the food area of the store, with a maximum flow rate of 30 m 3 day −1 (two batch phases every 12 h) and total nitrogen (TN) load of 135 mg L −1 . The IKEA ® store and the Department of Agriculture, Food, and Environment of the University of Catania have been in charge of running and monitoring the HSTW, respectively, from the time of the installation. Because of the high fluctuations of the flow and TN [23], the primary treatment was integrated with a hybrid TW in 2014. The system includes three in-series connected beds: HSTW, which allows for a reduction of organic matter and suspended solid (SS) concentrations, followed by two vertical subsurface flow TWs (V 1 -TW, V 2 -TW), which are designed to remove the wastewater organic matter, SS, and to further nitrify ammonia to nitrate. The HSTW bed has a surface area of about 400 m 2 (12 × 34 m), and is filled with 0.8-1.5 × 10 −2 m of volcanic gravel to a depth of 0.60 m.
The porosity and K s values of the original clean gravel were 0.41 and 19,466 m day −1 , respectively [18]. The HSTW, planted with Phragmites australis, is fed discontinuously with a daily effluent from the SBR plus effluent from the screening unit, which does not include the SBR, when the wastewater produced in the IKEA ® store exceeds a flow rate of 2.5 m 3 h −1 . Based on the total wastewater produced by IKEA, the annual volume of untreated wastewater feeding the HSTW was, on average, around 30% of the total volume, while during weekends and holidays, the daily volume of untreated wastewater (WW) was up to 50% of the total daily volume amount. Due to the high volumes of WW produced at the IKEA ® store at the end of 2016, the SBR was often bypassed, and a high organic load entered the hybrid TW system. Thus, it was already necessary to remove and replace the filter gravel close to the inlet area after three years of operation (in April 2017). A CR-1000 automatic weather station (Campbell Scientific, Logan, UT, USA) was installed close to the experimental plant to measure the temperature, wind speed and direction, rainfall, global radiation, and relative humidity. The vegetation harvest was carried out every year (at the end of January or at the beginning of February). During the investigation campaign, the vegetation groundcover was about 85%. In particular, the plant density was low in the area close to the inlet and, in some regions, randomly distributed along the middle path. For more details on the hybrid TW, see [23].

Ks Measurements in the Full-Scale HSTW Bed
The falling-head method was applied to determine the Ks values in the HSTW bed in February 2019. In particular, four falling-head infiltration tests were performed around each of the nine piezometers located in the bed (Figure 2) at the same depth (see Table 1). Two open-ended tubes, one impervious (IMP) and one pervious (P), and the corresponding equations adapted from [17] to various applicability conditions (Table 1) were used. In particular, the IMP permeameter (Table 1, Scheme 1, also known as the Standpipe) allowed us to evaluate the vertical Ks values using Equation (1). For isotropic gravels (clean or not clogged), this would not be a problem. However, given that the clogging process is likely not to be isotropic, the vertical Ks values obtained using this method could be different from the horizontal Ks values. Instead, the P permeameter (Table 1, Scheme 2) allowed us to evaluate both vertical and horizontal Ks, using Equation (2), as in [18]. Both schemes did not require drilling a borehole, as the pipe was merely pushed directly into the substrate. This saved time, compared to other schemes, and also minimized disturbance to the porous medium. Equation (1) is the following:

K s Measurements in the Full-Scale HSTW Bed
The falling-head method was applied to determine the K s values in the HSTW bed in February 2019. In particular, four falling-head infiltration tests were performed around each of the nine piezometers located in the bed (Figure 2) at the same depth (see Table 1). Two open-ended tubes, one impervious (IMP) and one pervious (P), and the corresponding equations adapted from [17] to various applicability conditions (Table 1) were used. In particular, the IMP permeameter (Table 1, Scheme 1, also known as the Standpipe) allowed us to evaluate the vertical K s values using Equation (1). For isotropic gravels (clean or not clogged), this would not be a problem. However, given that the clogging process is likely not to be isotropic, the vertical K s values obtained using this method could be different from the horizontal K s values. Instead, the P permeameter (Table 1, Scheme 2) allowed us to evaluate both vertical and horizontal K s , using Equation (2), as in [18]. Both schemes did not require drilling a borehole, as the pipe was merely pushed directly into the substrate. This saved time, compared to other schemes, and also minimized disturbance to the porous medium. Equation (1) is the following: where R and L are the radius and the submerged length (m) of the tube, respectively, as defined in [17]; and H 1 and H 2 are the water levels (m) in the permeameter cell corresponding to time t 1 and t 2 (s), respectively. Equation (2) is the following: where R mod and L mod are the radius (m) and the submerged length (m), respectively, as calibrated in [18]; and H 1 and H 2 are the water levels (m) in the permeameter cell corresponding to time t 1 and t 2 (s), respectively.
Water 2020, 12, x FOR PEER REVIEW  5 of 18 where R and L are the radius and the submerged length (m) of the tube, respectively, as defined in [17]; and H1 and H2 are the water levels (m) in the permeameter cell corresponding to time t1 and t2 (s), respectively. Equation (2) is the following: where Rmod and Lmod are the radius (m) and the submerged length (m), respectively, as calibrated in [18]; and H1 and H2 are the water levels (m) in the permeameter cell corresponding to time t1 and t2 (s), respectively. A small hole was dug in the granular medium to reach the water table, and then the permeameter was inserted using a mallet. A plastic water reservoir (6.6 L volume) with measurement units was assembled together with a ball valve to add water in a single-pulse mode, as required. A pressure probe (Sensor Technik Sirnach (STS), AG, Sirnach, Switzerland), connected to a laptop by means of a CR200-R (Campbell Scientific) data logger, was used to measure the variation of the water levels (H) within the measurement unit. A driver was added to the device to allow for the insertion of the pressure probe inside the steel permeameter ( Table 1). The pressure value at atmospheric pressure was checked before each measurement started. Four water level data per second were recorded for a duration of 30 s. The decrease in the water height inside the permeameter was monitored until the water reached the static water table. The best fit between the simulated and measured water levels was obtained by summing and minimizing the squared differences between the theoretical curve and that obtained in the field, following Equation (3): where Hobs is the height of the water table level measured inside the permeameter at time t during the test (m); and Hsim is the corresponding modeled data calculated by Equation (1) or Equation (2) (see Table 1), depending on the scheme used. This condition allowed us to estimate the value of Ks using an iterative, nonlinear procedure that makes use of the Excel Solver spreadsheet plug-in (Frontline Systems, Incline Village, NV, USA). A small hole was dug in the granular medium to reach the water table, and then the permeameter was inserted using a mallet. A plastic water reservoir (6.6 L volume) with measurement units was assembled together with a ball valve to add water in a single-pulse mode, as required. A pressure probe (Sensor Technik Sirnach (STS), AG, Sirnach, Switzerland), connected to a laptop by means of a CR200-R (Campbell Scientific) data logger, was used to measure the variation of the water levels (H) within the measurement unit. A driver was added to the device to allow for the insertion of the pressure probe inside the steel permeameter ( Table 1). The pressure value at atmospheric pressure was checked before each measurement started. Four water level data per second were recorded for a duration of 30 s. The decrease in the water height inside the permeameter was monitored until the water reached the static water table. The best fit between the simulated and measured water levels was obtained by summing and minimizing the squared differences between the theoretical curve and that obtained in the field, following Equation (3): where H obs is the height of the water table level measured inside the permeameter at time t during the test (m); and H sim is the corresponding modeled data calculated by Equation (1) or Equation (2) (see Table 1), depending on the scheme used. This condition allowed us to estimate the value of K s using an iterative, nonlinear procedure that makes use of the Excel Solver spreadsheet plug-in (Frontline Systems, Incline Village, NV, USA). Table 1. Schemes and equations tested for the K s measurements by the falling-head method, adapted from [17].

Permeameter Cell Scheme Symbols
Impervious, with the material inside (Standpipe) Applicability conditions: K s represents the vertical hydraulic conductivity of the soil Table 1. Schemes and equations tested for the Ks measurements by the falling-head method, adapted from [17].

Permeameter Cell Scheme Symbols
Impervious, with the material inside (Standpipe) Applicability conditions: Ks represents the vertical hydraulic conductivity of the soil R (0.05 m) and L (0.32 m) are the radius and the submerged length (m) of the tube, respectively, as defined in [17]. H1 and H2 are the water levels (m) in the permeameter cell corresponding to time t1 and t2 (s), respectively.
Pervious (for a total area of 3.11 × 10 −4 m 2 ), with the material inside. Applicability conditions: Ks represents both the horizontal and vertical hydraulic conductivity of the soil Rmod (0.049 m) and Lmod (0.20 m) are the radius (m) and the submerged length (m), respectively, as calibrated in [18]. H1 and H2 are the water levels (m) in the permeameter cell corresponding to time t1 and t2 (s), respectively. Lh (0.25 m) is the perforated permeameter length

Drainage Experiment
The drainage equation proposed in [20] and applied to HSTWs in [12] was used to find a synthetic Ks for the whole HSTW bed. The equation has the following form: where Pervious (for a total area of 3.11 × 10 −4 m 2 ), with the material inside. Applicability conditions: K s represents both the horizontal and vertical hydraulic conductivity of the soil Table 1. Schemes and equations tested for the Ks measurements by the falling-head method, adapted from [17].

Permeameter Cell Scheme Symbols
Impervious, with the material inside (Standpipe) Applicability conditions: Ks represents the vertical hydraulic conductivity of the soil R (0.05 m) and L (0.32 m) are the radius and the submerged length (m) of the tube, respectively, as defined in [17]. H1 and H2 are the water levels (m) in the permeameter cell corresponding to time t1 and t2 (s), respectively.
Pervious (for a total area of 3.11 × 10 −4 m 2 ), with the material inside. Applicability conditions: Ks represents both the horizontal and vertical hydraulic conductivity of the soil Rmod (0.049 m) and Lmod (0.20 m) are the radius (m) and the submerged length (m), respectively, as calibrated in [18]. H1 and H2 are the water levels (m) in the permeameter cell corresponding to time t1 and t2 (s), respectively. Lh (0.25 m) is the perforated permeameter length

Drainage Experiment
The drainage equation proposed in [20] and applied to HSTWs in [12] was used to find a synthetic Ks for the whole HSTW bed. The equation has the following form: where

Drainage Experiment
The drainage equation proposed in [20] and applied to HSTWs in [12] was used to find a synthetic K s for the whole HSTW bed. The equation has the following form: where K s (m min −1 ) is the only unknown term in the above equation, when applied to the drainage experiments. Measured cumulative drainage volume data were obtained first by closing the inflow and outflow valves to produce a horizontal water table and then by opening the outlet valve and collecting the drainage volume in a tank. The inflow valve remained closed for the whole duration of the test. The variation of the water levels within the tank was measured by the pressure probe positioned in the tank and connected to a laptop, as described above. The water level data, recorded inside the tank every 5 min, were converted into drainage volumes. The drainage experiment was carried out once in February 2019. A range of conductivity values were applied to Equation (4), until the resulting outflow curve matched the measured outflow. The value of K s that produced the closest fit and gave the simulated cumulated outflow closest to the measured one was taken as the hydraulic conductivity of the substrate for the particular drainage test.

Tracer Tests
A tracer test was conducted in February 2019 by introducing an impulse of NaCl into the HSTW inlet at time zero. NaCl was chosen for its efficiency and low-cost characteristics, but also because it is not organic, and it was already successfully used in the same system [23]. The amount of tracer (32 kg 100 L −1 for a total volume of 200 L and a pulse duration of about 10 min) was chosen in order to reach an average benchmark concentration at least ten times that of the background concentration, as suggested in [26]. The measured EC (µS cm −1 ) values, after subtraction of the background value (between 1600 and 1800 µS cm −1 ), were then converted into NaCl concentrations (mg L −1 ) by a linear calibration curve (R 2 = 0.99). The pumping discharge of the solution was the same as the normal water flow in the HSTW system. The solution was prepared in a bucket, in which the tracer was added and mixed until a uniform concentration was reached. A pulse-inject of the tracer solution was then added with a pump into the inlet zone of the HSTW system. The fluid electric conductivity (EC) was measured by ten conductivity probes (delta OHM-HD 2106.2, DeltaOhm, Padova, Italy). One of them was located at the outlet, and the others were located inside the nine piezometers in the HSTW bed. The EC data were recorded by each probe every 15 min. The HSTW outflow volumes were measured using a flow measurement device for the whole test. No rain occurred during the test.
The NaCl injection into the HSTW inlet provided information about the efficiency and detention times in the system. In fact, in ideal flow patterns, the water particles move at the same velocity, reaching the outlet together. In this case, a tracer impulse would also exit as an impulse (a sharp spike of concentration). It is clear from numerous studies that TWs are neither plug flow nor well mixed, and due to preferential flow channels, vertical stratification occurs in gravel beds, with more significant flows arising at lower levels in the system [27][28][29]. Thus, also in the present case, the response to the impulse tracer input was a time-delayed, bell-shaped curve, called the retention time distribution (RTD). Thus, the analysis of this curve allowed for the derivation of critical parameters characterizing the hydraulic behavior of the system.
First, it was necessary to evaluate the validity of the test, checking that the tracer was recovered nearly in its entirety at the wetland outlet. To achieve this, the relative tracer mass recovery (%) was calculated as M out M in × 100 (5) where: M out is the total recovered tracer mass at the outlet (g); M in is the total tracer introduced into HSTW (g); and M out was determined in accordance with [30]: where Q(t) is the outflow rate at time t (m 3 h −1 ); C(t) is the tracer concentration at time t at the outlet (mg L −1 ); and t is the sample time (h).
Then, the retention time distribution (RTD), which represents the various time fractions of the water spent in the reactor and hence the contact time distribution for the system was analyzed.
In general, the RTD is the probability density function for the residence times in a wetland. For an impulse input of a tracer into a steadily flowing system, the time function is defined by The first numerator is the mass flow of the tracer in the wetland effluent at any time t, after the time of the impulse addition. The first denominator is the sum of all the tracers collected and should thus equal the total mass of the injected tracer.
Thus, the mean tracer detention time (t m ) is presumed to be the actual mean detention time and can be calculated as the first moment of the curve (the centroid of the C-curve) as where t is the sample time that, after the definition of the probability density function, can be expressed as A wetland may have internal excluded zones that do not interact with the flow [4], such as the volume occupied by biomass (roots, rhizomes, etc.). In a steady-state system without excluded zones, the tracer detention time (t m ) equals the nominal residence time (t n ), which is defined as where V (m 3 ) is the volume of the water in the system; and Q (m 3 h −1 ) is the water flow rate through the system.
Following [31], a parameter for evaluating the hydraulic efficiency of the system can be simply defined as where t p is the time that the maximum concentration of tracer occurred (h). This measure has the advantage of being readily derived from the RTD, and it does not have the problem related to the calculation of t m [31]. The hydraulic efficiency can be considered good if λ > 0.75; it is satisfactory if 0.5 < λ < 0.75; and it is poor if λ ≤ 0.50. A second parameter that can be determined directly from the residence time distribution is the variance (σ 2 , Equation (12)) (or the square of the standard deviation): where σ 2 is the second moment of the curve and characterizes the spread of the tracer response curve concerning the mean of the distribution (t m ).

Wastewater Quality Characterization
Wastewater quality was monitored one/two times per month, from April 2016 to February 2019. The sampling points were located at the hybrid TW inlet and at the three outlets of each TW bed (HSTW, V 1 -TW, and V 2 -TW). Standard methods [32] were applied for wastewater quality analysis, which included total suspended solids (TSS, mg L −1 ) at 105 • C, biological oxygen demand after five days (BOD 5 , mg L −1 ), chemical oxygen demand (COD, mg L −1 ), phosphate (PO 4 , mg L −1 ), ammonia nitrogen (N-NH 3 , mg L −1 ), total nitrogen (TN, mg L −1 ), nitrate nitrogen (N-NO 3 , mg L −1 ), and Escherichia coli (E. coli, CFU 100 mL −1 ). In order to characterize the HSTW influent, the following parameters were calculated [2].
The maximum areal organic loading rate (OLR) (g BOD 5 m −2 day −1 ) is defined as The hydraulic loading rate (HLR) (mm day −1 ) is defined as The total suspended solids rate (TSSLR) (g TSS m −2 day −1 ) is defined as where Q (m 3 day −1 ) is the water inflow rate through the system; TSS (g day −1 ) is the amount of TSS, based on the water inflow rate through the system; BOD 5 (g day −1 ) is the amount of BOD 5 based on the water inflow rate through the system; and A (m 2 ) is the wetland surface area. Table 2 and Figure 3 report the K s mean values (indicated as K s values below) around each of the nine piezometers inside the HSTW bed and their standard deviation (SD) observed in 2019 and in the previous campaign [18]. The K s values measured in piezometer 1, 2, and 3, located in the first transect close to the inlet, were the lowest. A very low K s value was also observed in piezometer 8 (Tables 2  and 3). Consequently, in 2019 the lowest values of K s were observed in the inlet transect (660 m day −1 ) ( Table 3) and along the middle path (3682 m day −1 ) ( Table 4). The K s values found along the transects located in the middle part of the bed and close to the outlet were similar (Table 3).

Drainage Experiment
After the opening of the outlet valve, the outflow of the HSTW bed measurements were recorded every 5 min and ended after 1250 min. The K s value that gave the best fit between the observed and simulated cumulative outflows was 3880 m day −1 (Figure 4). The difference between the observed and simulated cumulative flows was 0.1 m 3 . Figure 4 shows the observed and simulated cumulative outflow curves during the drainage experiment.

Tracer Tests
The NaCl concentrations obtained at the outlet of the HSTW bed during the tracer test experiment are shown in Figure 5. The overall percentage of tracer recovery was 94.5% (M out = 60 × 10 3 g), after a duration of 195 h. The high percentage of tracer recovery indicated that NaCl was marginally affected by the clogging matter. The mean residence time calculated using the moment analysis method was 88 h (with σ 2 = 46 h). During the tracer test, the average flow rate was 1.

Tracer Tests
The NaCl concentrations obtained at the outlet of the HSTW bed during the tracer test experiment are shown in Figure 5. The overall percentage of tracer recovery was 94.5% (Mout = 60 × 10 3 g), after a duration of 195 h. The high percentage of tracer recovery indicated that NaCl was marginally affected by the clogging matter. The mean residence time calculated using the moment analysis method was 88 h (with σ 2 = 46 h). During the tracer test, the average flow rate was 1.       (Table 5) confirmed that a large number of customers visit the IKEA ® store during weekends and holidays. In particular, on these busy days, the hybrid TW often received wastewaters directly from the screening unit (up by 40-50% from the total daily HSTW influent) and therefore had a lower quality. The hybrid TW units provided a very high mean efficient reduction of TSS, COD, and BOD5 especially in the HSTW, thus reducing the clogging problem in V1-TWs and V2-TWs and allowing for the limits fixed by the Italian law to be respected. The TN reduction was very high, confirming that both processes (nitrification and denitrification) were efficient.     (Table 5) confirmed that a large number of customers visit the IKEA ® store during weekends and holidays. In particular, on these busy days, the hybrid TW often received wastewaters directly from the screening unit (up by 40-50% from the total daily HSTW influent) and therefore had a lower quality. The hybrid TW units provided a very high mean efficient reduction of TSS, COD, and BOD5 especially in the HSTW, thus reducing the clogging problem in V1-TWs and V2-TWs and allowing for the limits fixed by the Italian law to be respected. The TN reduction was very high, confirming that both processes (nitrification and denitrification) were efficient.  Table 5 shows the mean values of the physical-chemical parameters (and their standard deviations) evaluated by 26-38 samples for each parameter in different sections of the hybrid TW system from April 2016 to February 2019. The limits imposed by Italian regulations for wastewater discharge into water bodies are also shown (Table 3, annex 5, III part of [33]). A high variability of pollutant concentrations at the inlet of the hybrid TW evaluated by 26-38 samples for each parameter collected during the period from April 2016 to February 2019 (Table 5) confirmed that a large number of customers visit the IKEA ® store during weekends and holidays. In particular, on these busy days, the hybrid TW often received wastewaters directly from the screening unit (up by 40-50% from the total daily HSTW influent) and therefore had a lower quality. The hybrid TW units provided a very high mean efficient reduction of TSS, COD, and BOD 5 especially in the HSTW, thus reducing the clogging problem in V 1 -TWs and V 2 -TWs and allowing for the limits fixed by the Italian law to be respected. The TN reduction was very high, confirming that both processes (nitrification and denitrification) were efficient.

Wastewater Quality Characterization
Moreover, because of the effluent recirculation to primary treatment, the nitrate at the hybrid TW outlet was lower than the Italian discharge threshold. Instead, the mean TP concentration was higher than the limit, which is probably due to the filter medium composition. A very significant mean reduction of E. coli was found at the hybrid TW (3 log unit). The mean annual values (and their standard deviations) for OLR, HLR, and TSSLR evaluated at the inlet of the HSTW are reported in Table 6. The mean annual water flow rate varied from 24 to 32 m 3 day −1 , with TSSLR, HLR, and OLR values in the range of those suggested for correct HSTW design [2]. TSSLR and HLR increased from 2016 to 2018.

Discussion
To better understand the hydraulic behavior of the HSTW, K s measurements were conducted in February 2019 using the newly implemented scheme (a pervious permeameter and calibrated equation) and were compared with data obtained in April 2018 [18], and with data obtained for clean gravel (K s = 19,466 m day −1 ) ( Table 2). It is important to highlight that the wastewater characteristics (flow and quality), vegetation, and the main features of the K s experimental campaigns of 2018 and 2019 were very similar. As expected, the spatial evolution of the clogging since the beginning of the operation period (2014) and during the observation period (2018-2019) was not uniform within the HSTW bed. It was more severe in the area close to the inlet [15,34,35] and along the central path. Considering the temporal evolution of the phenomenon, very high reductions of K s values (up to 89%) were already observed in 2018 after four years of operation for clean gravel. These reductions were generally more robust in 2019 (up to 99%) for the nine piezometers (Table 2). Only the K s value for piezometer 9 did not increase in 2019 (K s = 7285 m day −1 ). Regardless, the variation in the year 2018 (976 m day −1 ) was lower than that of the SD observed in 2019 (±1223 m day −1 ). The lowest K s value found in the inlet transect was also in 2018 (3221 m day −1 ), but a very high reduction occurred in 2019 ( Table 3). The K s values along the central path 2-5-8 observed in 2019 were lower than those noted in 2018 ( Table 4). The K s values that occurred along the lateral paths in 2018 and 2019 were similar ( Table 4). The slight K s reduction found in the outlet zone in 2018 (also observed in [36]) was also confirmed in 2019. The values of K s that varied from a minimum of 225.8 m day −1 (±125.2) to a maximum of 7468.3 m day −1 (±517.1) at piezometer 2 and 6, respectively, in 2019, were in the range of those found by other authors for HSTW gravels of a similar size and porosity [15,19,37]. Some studies [6,7] reported lower K s values (maximum values at the outlet area: 810 m day −1 ) using the falling-head method with a pervious permeameter. Moreover, recent studies confirm the difficulty of measuring the hydraulic conductivity in HSTWs [15,22,38], especially when the clogging is in the beginning phase and the growing vegetation makes the substrate nonisotropic.
Notwithstanding the information acquired from punctual K s measurements at nine points inside the HSTW, it was difficult to evaluate the K s value of the entire system, which cannot simply be the arithmetic mean, which, in this case, is 4425 m day −1 . Thus, for this evaluation, the drainage method proposed in [20] was used. The K s value corresponding to 3880 m day −1 found for the HSTW bed was similar to the arithmetic mean value of the K s measurements obtained inside the bed. The K s value was obtained after a calibration of Equation (4). The HSTW was put out of service for the duration of the test (1250 min in this case). A K s value of 3456 m day −1 was found in [12] by applying the drainage equation method in a pea gravel bed with a thickness of 0.5 × 10 −2 m, after a 2-year operation period. The drainage method also allowed for the evaluation of the porosity reduction of the gravel. In particular, the observed cumulative outflow revealed a porosity variation of the HSTW gravel (with a size in the range from 0.8 × 10 −2 to 1.5 × 10 −2 m) from 0.41 to 0.28, which supports the reduction of the available water volume in the bed due to the clogging. A similar porosity reduction in HSTW beds was observed in [12]. In particular, a variation from 0.33 to 0.27 and 0.36 to 0.33 was found for pea gravel (with a size of 0.5 × 10 −2 m) and coarse gravel (with a size in the range of 2 × 10 −2 to 4 × 10 −2 m), respectively, after a 2-year operation period. Note: TSS, total suspended solids (mg L −1 ); BOD 5 , biological oxygen demand after five days (mg L −1 ); COD, chemical oxygen demand (mg L −1 ); PO 4 , phosphate (mg L −1 ); TN, total nitrogen (mg L −1 ); N-NH 3 , ammonia nitrogen (mg L −1 ); N-NO 3 , nitrate nitrogen (mg L −1 ); E. coli (CFU 100 mL-1). a [33]; b Limit for discharge into surface water bodies. As expected, the spatial evolution of clogging since the beginning of the operation period (2014) and during the observation period (2018-2019) was not uniform inside the HSTW bed, and was more severe in the area close to the inlet [15,34,35]) and along the central path. Considering the temporal evolution of the phenomenon, very high reductions of the K s values were observed already after four years of operation in 2018 (up to 89%) for clean gravel. These reductions were generally more robust in 2019 (up to 99%) for the nine piezometers (Table 2). Only the K s value for piezometer 9 did not increase in 2019 (K s = 7285 m day −1 ). Regardless, the variation for the year 2018 (976 m day −1 ) was lower than that of the SD observed in 2019 (±1223 m day −1 ). The lowest K s value found in the inlet transect was also in 2018 (3221 m day −1 ), but a very high reduction occurred in 2019 ( Table 3). The K s values along the central path 2-5-8 observed in 2019 were lower than those noted in 2018 ( Table 4). The K s values that occurred along the lateral paths in 2018 and in 2019 were similar ( Table 4). The slight K s reduction found in the outlet zone during 2018 (also observed by [36]) was confirmed in 2019. The values of K s that varied from a minimum of 225.8 (±125.2) to a maximum of 7468.3 m day −1 (±517.1) at piezometer 2 and 6, respectively, in 2019 were in the range of those found in other research for HSTW gravels of a similar size and porosity [15,19,37]. Some studies [6,7] reported lower K s values (maximum values at the outlet area: 810 m day −1 ) using the falling-head method with a pervious permeameter. Moreover, recent studies confirm the difficulty of measuring the hydraulic conductivity in HSTWs [15,22,38], especially when clogging is in the beginning phase and the growing vegetation makes the substrate nonisotropic.
Several authors reported that a comparison of tracer tests performed during different periods is essential for understanding the spatial and temporal evolutions of clogging. Therefore, the results of the tracer test carried out in 2019 were analyzed in relation to those acquired for the same HSTW system in the past [23], along with the results of the K s and drainage experiments carried out in 2019. The integrity of the tracer test results was guaranteed by the duration of the test, which lasted for 195 h in 2019 and is twice the nominal residence time [39]. The mean time was higher than the nominal residence time of about nine hours, thus confirming the presence of dead or stagnant zones [40] and indicating that the flow may not include the entire volume of the HSTW bed [4,41]. In any case, the hydraulic efficiency observed in 2019 was good (λ = 0.80), following [31]. The passage of the tracer monitored at the outlet revealed a uniform RTD, with a single peak of the NaCl concentration of about 480 mg L −1 occurring after 62.25 h from the beginning of the test, probably indicating a preferential path through the wetland [19]. Based on the K s measurement results, the preferential path should be the central one (2-5-8), because the highest reductions of K s were measured during the observation period. The preferential path inside the HSTW bed seems to be also confirmed by the results of the tracer tests carried out in [23]. The EC (mS cm −1 ) values were higher along the central path than along the lateral paths during both tracer tests carried out in February and May 2017 [23]. Moreover, the EC maximum values observed at piezometers 2 and 5 in the present study (6950 and 2545 mS cm −1 , respectively) were lower (and the peak time was delayed) than those achieved in [23] before the restoration measurement in the inlet zone and higher (and the peak time was accelerated) than those observed after the restoration measurement. These differences can be explained by the fact that the clogging in the HSTW bed is currently in an intermediate situation, in comparison to those reported in [23], indicating that it will soon require a restoration measure (i.e., removal and replacement of the filter gravel close to the inlet area). Thus, our findings support the idea that a comparison with previous tracer tests in the same system is fundamental for correctly interpreting spatial and, obviously, temporal hydraulic changes in HSTWs, as emphasized by several authors [19]. In [21], it was found that in an HSTW bed, an area which exhibited a substantial preferential flow path had become a large dead zone after a 2-year operation period.
Another important finding of the present study is that K s values and the flow distribution obtained using widely applied methods, such us the falling-head method (with a newly implemented permeameter) and traditional tracer tests, are consistent with each other and with the K s values obtained by using the drainage equation method proposed in [20]. The sensitivity in assessing the clogging phenomenon inside wetlands using different K s measurements inside the bed and path flows (obtained by tracer tests) was already highlighted in several studies [19,23]. The use of a drainage equation method that can provide an estimate of the hydraulic behavior of an entire system by determining a single value of K s can further help in the analysis of the complex hydrology of HSTW beds. Unfortunately, in the latter case, it is necessary to put the HSTW out of service. This limits the feasibility of the test for large HSTWs.
Moreover, considering the qualitative aspects of wastewater, the hybrid TW (SBR and TW) treating up to 50 m 3 day −1 water flow and characterized by the TSSLR, HLR, and OLR design parameters, as suggested in the literature on HSTWs, was able to function within Italian law standards, with only partial clogging of the HSTW bed.

Conclusions
Bed clogging of horizontal subsurface treatment wetlands (HSTWs) is a complex phenomenon that must be continuously monitored, especially in systems that have already been affected with clogging during their lifetime, in order to apply restorative measurements and prevent negative effects on treatment efficiency. In the HSTW of the hybrid TW system functioning as a secondary wastewater treatment system of the IKEA ® store, it was already necessary to remove and replace the filter gravel close to the inlet area after three years of operation (in April 2017). The application of traditional and innovative methods to assess the clogging allowed for the quantification of the phenomenon in the HSTW, as well as the identification of the most time-efficient and easily repeatable way to quantify the phenomenon. First, it must be noted that the K s measurement (obtained by a newly implemented scheme), tracer test, and drainage equation results were consistent with each other in highlighting a partial clogging in the inlet area of the HSTW bed of the hybrid TW system after 5 years of operation, despite the fact that the filter gravel was replaced in 2017. The K s values were very low (<1000 m day −1 ) in the transect close to the inlet and along the central flow path. The drainage method revealed a high reduction of porosity in the HSTW gravel during the operation period. The RTD was uniform, with a single peak of the NaCl concentration, and the mean residence time was higher than the nominal time. This confirmed that there was nonuniform flow through the HSTW and the presence of stagnant zones. Moreover, the comparison of the results with tests carried out previously for the same system suggests the need to perform these measurements at least once a year and to plan HSTW restorative measures soon.
Moreover, considering the main advantages and disadvantages of the methods applied, it must be noted that the newly implemented scheme (with the pervious permeameter) reduces the number of K s measurement points along the vertical dimension, and only measurements along the paths and transects are needed to characterize K s . On the other hand, the drainage equation method, allowing for the identification of one K s value for the whole bed at once, requires the system to be out of service for the entire duration of the test. Finally, tracer tests allow for an understanding of the effect of the clog matter on the flow through a porous medium, rather than assessing the clogging severity. In conclusion, the results highlighted that, especially for large systems, the combination of the newly implemented scheme with traditional tracer tests is the most efficient approach for more thoroughly understanding clogging development in HSTWs.
Despite the partial clogging of HSTW, the hybrid TW was able to meet the Italian law standards for wastewater discharge into water bodies.