Treatment of Agricultural Drainage Water by Surface-Flow Wetlands Paired with Woodchip Bioreactors

: Nutrient losses from agricultural ﬁelds have long been a matter of concern worldwide due to the ecological disturbance this can cause to surface waters downstream. In this paper a new design concept, which pairs a surface-ﬂow constructed wetland (SFW) with a woodchip bioreactor (WB), was tested in relation to its capacity to reduce both nitrogen (N) and phosphorus (P) loads from agricultural tile drainage water. A nutrient mass balance and a comparative analysis were carried out together with statistical regressions in order to evaluate the performance of four SFW + WBs under di ﬀ erent catchment conditions. We found marked variations between the systems in regard to hydraulic loading rate (0.0 to 5.0 m / day) and hydraulic retention time (1 to 87 days). The paired system worked as nutrient sinks throughout the study period. Total N and total P removal e ﬃ ciencies varied from 8% to 51% and from 0% to 80%, respectively. The results support the use of the new design concept for nutrient removal from tile-drained agricultural catchments in Denmark as part of national management plans, with the added advantage that smaller areas are needed for construction (0.1% to 0.2% of the catchment area) in comparison to standalone and currently used SCWs (~1%).


Introduction
Most agricultural areas represent a major source of nitrogen (N) and phosphorus (P) to surface waters [1][2][3]. Excessive nutrient loads to marine and freshwater ecosystems may cause eutrophication [4,5] followed by the loss of biodiversity [6,7]. In this context, the EU Water Framework Directive 2000/60/EC requires that nutrient concentrations do not exceed levels that disturb the normal functioning of the natural ecosystem thus putting emphasis on the need for substantial removal. In Denmark, agriculture covers 61% of the country area (Statistics Denmark, Statistical yearbook 2017), in which approximately 50% is tile-drained [8]. The N and P losses via tile drainage amount to 45-60% (~22,000 t N/year) and approximately 33% (~400 t P/year) of the total losses [9]. This represents a substantial hazard to surface waters downstream, as natural removal processes in the landscape are bypassed by the tile drainage networks, which function as a conduit for nutrients. Therefore, current national efforts are searching for effective measures at the landscape scale. These can intercept the drainage networks and reduce the nutrient load to surface waters, in order to prevent eutrophication and meet the environmental objectives of the EU Water Framework Directive.
Surface-flow constructed wetlands (SFWs) are deemed cost-effective in the removal of nutrients from agricultural drainage discharge [10], and are often reported as nutrient sinks [11][12][13][14][15][16][17][18][19][20]. However, some studies reported SFWs as P sources when receiving P predominantly in dissolved form [20,21]. This is generally attributed to removal processes for nitrate being superior to those for phosphate [22]. In these systems, removal of N and P loads is dependent on the nutrient load and hydraulic retention time (HRT), thus being variable when receiving event-driven drainage discharge [23][24][25]. Nutrient removal rate commonly shows a strong positive correlation to the nutrient load, while effective nutrient removal (%) largely depends on sufficiently long HRT [26][27][28]. Removal of N mainly occurs by denitrification in the anaerobic sediments and is, to a high extent, controlled by carbon availability and temperature [24,29]. Retention of P, in turn, mainly occurs by sedimentation of particle-bound P (PP) and sorption of dissolved reactive forms (e.g., orthophosphate) to reactive sites in the soil/sediments [30]. Moreover, micro and macro-organisms (e.g., plants) can also reduce the nutrient load by biological uptake. However, storage of organic compounds is generally a short-term process leading to nutrient release during mineralization [29,30]. According to Werner and Kadlec [31], the hydraulic loading rate (HLR) also plays a major role in controlling P retention and influences the water flow patterns, which regulate the HRT distribution across the system. However, previous studies failed to show a strong correlation between TP retention efficiency and HLR [26,32].
Heavy precipitation events lead to high HLR and short HRT. This can hinder denitrification and P sorption, subsequently limiting the performance of an SFW. Attempts to enhance the nutrient removal processes in this situation may include the paring of an SFW with a woodchip bioreactor (WB). Removal of N from agricultural drainage discharge has been widely studied in the latter systems, where performance has been reported as markedly higher [33][34][35][36][37]. Moreover, some studies even documented higher N removal rates in WBs than in SFWs [38,39]. These results can be attributed to the anaerobic medium and large carbon availability in WBs, which favors denitrification [40]. These systems can ensure N removal for long periods (e.g., for more than five years) as woodchip degradation is rather slow under anaerobic conditions, with an estimated half-life between 5 and 37 years and likely depending on the redox conditions and denitrifying activity in the system [35]. Among other factors, the woodchip degradability and water temperature largely regulate the N removal rate in WBs [40,41]. In addition, increasing the inlet N concentration and HLR may enhance the removal rate, providing that there is carbon availability and proper temperature. Effective N removal (%), on the other hand, is reported to increase under longer HRTs [42]. Despite the scarcity of studies on P retention by WB, Gottschall et al. [43] demonstrated that these systems not only remove N, but also retain P from agricultural drainage discharge. Moreover, Choudhury et al. [44] demonstrated that P retention was significant in effluents containing mostly PP (83% of total P), while the retention of dissolved P was minimal.
In Denmark, Mendes et al. [18] reported SFWs as P sinks, retaining from 0.3 to 10.5 g/m 2 /year, resulting in 41% to 51% P retention over a three-year period. Bruun et al. [45], in turn, reported N removal in Danish WBs with horizontal or vertical flow hydraulic design from 8.0 to 10.5 g/m 2 /day (water flow = 1.8 L/s), resulting in 45-57% N removal. However, a design concept consisting of pairing an SFW with a WB, thus potentially enhancing the nutrient removal processes in a single system remains untested. The construction of the new design concept utilizes less area in relation to the catchment area (approximately 0.1% to 0.2%) in comparison to currently used SFWs in Denmark (~1%). Therefore, this study aimed to assess the annual and seasonal performance in the removal of N and P by a newly developed design concept, consisting of paired systems (SFW+WB) receiving agricultural tile drainage water in Denmark. Moreover, a performance comparison was carried out between a SFW+WB and a SFW at a common location, thus receiving identical nutrient concentrations at the inlet. We hypothesized that the SFW+WB would outperform in comparison with systems with similar or larger size in relation to N and P removal. The results presented herein will ultimately help authorities to evaluate the viability of the paired systems as part of management plans to reduce nutrient loads from Danish agricultural catchments to surface waters.

Figure 1.
Overview of the investigated surface-flow wetlands (SFWs) paired with woodchip bioreactors (WBs) at (a) Aakaer, (b) Ondrup, (c) Ryaa, and (d) Gjol, the latter also including the standalone SFW. The main components are indicated as sedimentation basin (SB), woodchip bioreactor (WB), clarification basin (CB), in addition to inlet (in) and outlet (out). The SFW+WBs have the following 1-m deep components in the order of the water flow direction: (i) a sedimentation basin (SB), (ii) a woodchip bioreactor (WB) consisting of raw woodchip material having a diameter of approximately 95 mm (Ø land Entreprenørforretning, Aabybro, Denmark), and (iii) a clarification basin (CB). Tile drainage water is collected from adjacent agricultural fields and pumped (E15A-2/3, KW PUMPER Wünche A/S, Vae rlose, Denmark) into the SB to allow sedimentation of suspended particles and particle-bound P (PP). The water then  The SFW+WBs have the following 1-m deep components in the order of the water flow direction: (i) a sedimentation basin (SB), (ii) a woodchip bioreactor (WB) consisting of raw woodchip material having a diameter of approximately 95 mm (Øland Entreprenørforretning, Aabybro, Denmark), and (iii) a clarification basin (CB). Tile drainage water is collected from adjacent agricultural fields and pumped (E15A-2/3, KW PUMPER Wünche A/S, Vaerlose, Denmark) into the SB to allow sedimentation of suspended particles and particle-bound P (PP). The water then infiltrates through the WB, representing the carbon source to activate denitrification. In relation to P, the WB acts as a medium for P immobilization and deposition of PP and iron in the matrix voids. The raw woodchip material is partly degraded with time and the remaining part is dismissed as waste. The CB follows the WB, and is intended to oxygenate the anoxic effluent leaving the WB. The areas and volumes of each component vary between the SFW+WBs ( Table 1). The area ratios of the paired systems to their catchments range from 0.11% to 0.19%. An SFW was also constructed in Gjol during 2015 and is located south of the SFW+WB (Figure 1d). The goal was to compare the nutrient removal performance of both systems with identical nutrient concentrations at the inlet. The area ratio of the SFW to its catchment is 0.49%.

Water Sampling and Chemical Analyses
Electromagnetic flowmeters (OPTIFLUX 2000 Sensor, Krohne, Duisburg, Germany) were installed at the inlet of the SFW+WBs, and allowed continuous recording of the inflow drainage discharge (Q). The SFW+WB in Ondrup was equipped from September 2017 with an additional flowmeter, which measured any discharge (Q overflow ) bypassing the system through a parallel pipe connecting the inlet with the outlet. The inlet flowmeter in Gjol is located in a common distribution well, from which drainage water subsequently splits into the SFW+WB and the SFW (water flow ratio~1:2.5). The HLR was the ratio between Q and the SFW+WB area, while the HRT was the ratio between the SFW+WB volume and Q. Automatic water samplers (3700 Sampler, Teledyne ISCO, Lincoln, NE, USA) were installed at the inlet and outlet of the SFW+WBs. However, the study site in Gjol, was equipped with automatic water samplers at the common distribution well (inlet), at the outlet of the SB and WB, and at the outlet of the SFW. Water samples of 30 mL were collected automatically every hour and pooled into a daily sample in a 1-L polyethylene bottle. Samples were collected every 21 days, transported to the lab, and stored at 2 • C for approximately 1 month before chemical analysis. Depending on the variations of the drainage discharge hydrograph, analysis was carried out either individually or as a composite sample. Water samples collected in periods of highly variable water flow were analyzed individually to measure nutrient concentrations during peak-flow events, while those collected in periods with fairly steady water flow were pooled and analyzed as composite samples, as the variation in nutrient concentration was expected to be lower [18]. Between the end of June and the start of September 2018, data were not recorded in Ryaa due to dry weather conditions and no tile drain discharge at the inlet. For total N (TN) analysis, unfiltered water samples were initially autoclaved for 30 min at 121 • C in an alkaline peroxydisulphate solution. The digested samples were subsequently acidified and their adsorbance was measured using an auto analyzer (AutoAnalyzer 3, SEAL, Norderstedt, Germany). Total P (TP) was determined in a spectrophotometer (Spectronic Hellios Alpha, Thermo ScientificTM, Dartford, UK) at 890 nm with ascorbic acid (5%), after 30 min acid persulfate (5%) digestion in autoclave at 121 • C [46,47]. Chemical analysis for other pollutants was not carried out as it is not of concern under the investigated conditions. Water sampling and chemical analysis followed previous established methodologies used for nutrient mass balance studies in Denmark [18,23,24].
The monitoring at Aakaer and Ondrup started in 2012 and 2010, respectively. However, some technical adjustments to the original design were made in 2016. Thus, the data monitored up to 2016 were excluded from the analyses carried out in this paper. Accordingly, a new phase of monitoring started in 2016 under the Future Cropping project (https://futurecropping.dk/). The monitoring at Ryaa, in turn, also started in 2016, although the system was constructed earlier. Since 2018, monitoring has also included the SFW+WB and SFW in Gjol. The term hydrological year in this paper refers to the period from the 1st of August to the 31st of July. Moreover, the term winter half year refers to the period spanning from October to the end of March, while summer half year refers to the remaining annual period (April to the end of September).

Statistical Analysis
Simple regression analysis was performed between (i) removal rate and load, (ii) removal efficiency and HLR, and (iii) removal efficiency and HRT, in which both TN and TP data were tested. The regressions were either linear or quadratic, depending on the explanatory power of the model. Moreover, multiple regressions were used to test the effect of both HLR and TN or TP concentration (independent variables) on (i) the removal rate and (ii) removal efficiency of TN and TP (dependent variables). The statistical analyses were carried out in Microsoft Excel 14.0 (Microsoft, Redmond, WA, USA) using monthly values (n = 35-36 in Aakaer, Ondrup and Ryaa; n = 11 in Gjol), in which HRT and the concentrations of TN and TP were averages. The tests used a 95% confidence level.

Hydrological Performance
The drainage discharge (Q) for the different hydrological years is summarized in Table 2. Aakaer received the highest annual Q, followed by Ondrup and Ryaa. The HLR variation was similar in the three SFW+WBs ( Figure 2). Aakaer (Figure 2a) received the highest loads, with peaks ranging from 3 to 5 m/day. Ondrup and Ryaa (Figure 2b,c) presented lower values, with peaks ranging from 1 to 3 m/day. Three periods of intense HLR were observed during the winter half year, while lower and steadier HLR was observed during the summer half year. The variation of the HRT is opposite to that of the HLR. A number of three high HRT periods can be identified, corresponding to the summer half year (low daily values of HLR). Average daily HRT during the summer half year was equal to 8, 7, and 87 days, while during the winter half year it was equal to 1, 2, and 5 days for Aakaer, Ondrup, and Ryaa, respectively. This relates to the differences in water volume in the paired systems, corresponding to 906, 1062 and 1208 m 3 in Aakaer, Ondrup and Ryaa, respectively.

Nutrient Loads and Removal
Daily TN concentrations ranged from 1 to 23 mg/L and from 0 to 16 mg/L at the inlet and outlet, respectively. Monthly discharge-weighted TN concentration varied similarly to HLR, where Aakaer ( Figure 3a) presented the highest values followed by Ondrup and Ryaa (Figure 3b,c). Seasonal variations in TN concentration also occurred, with higher and lower values during the winter and summer half year, respectively. Daily TP concentrations ranged from 0.01 to 1.00 mg/L and from 0.00 to 1.07 mg/L at the inlet and outlet, respectively. Monthly discharge-weighted TP concentration for Aakaer and Ondrup (Figure 3a,b) were comparable with clear peaks. However, these concentrations in Ryaa (Figure 3c) were lower and steadier, so that no clear seasonal variation was visible.

Nutrient Loads and Removal
Daily TN concentrations ranged from 1 to 23 mg/L and from 0 to 16 mg/L at the inlet and outlet, respectively. Monthly discharge-weighted TN concentration varied similarly to HLR, where Aakaer (Figure 3a) presented the highest values followed by Ondrup and Ryaa (Figure 3b,c). Seasonal variations in TN concentration also occurred, with higher and lower values during the winter and summer half year, respectively. Daily TP concentrations ranged from 0.01 to 1.00 mg/L and from 0.00 to 1.07 mg/L at the inlet and outlet, respectively. Monthly discharge-weighted TP concentration for Aakaer and Ondrup (Figure 3a,b) were comparable with clear peaks. However, these concentrations in Ryaa (Figure 3c) were lower and steadier, so that no clear seasonal variation was visible. The differences between inlet and outlet for monthly average discharged-weighted nutrient concentrations were generally higher during the summer half year and highest for the hydrological year 2018/2019 (Figure 3). For Aakaer, TN removal efficiency fluctuated between 9% and 16% during 2016-2019, being lower compared to previous hydrological years ( Table 2). Removal efficiency of TP here was, on the other hand, the highest compared to the other SFW+WBs, and reached 79% in 2018/2019. For Ondrup, TN removal efficiency varied between 16% and 32%, while for TP it varied between −23% and 17%. For Ryaa, TN removal efficiency varied between 8% and 32%, with minimal effect in 2016/2017 and maximum in 2018/2019. Removal efficiency of TP varied between 0% and 44%, and was generally high, despite the low inflow and outflow TP concentrations.
The simple regressions carried out in this study demonstrated that HLR correlated inversely to TN removal efficiency in the SFW+WBs (Figure 4). However, in Aakaer (Figure 4a) a quadratic model explained nearly half of the variation in TN removal efficiency and indicated that an increase may be achieved after a certain HLR threshold. When correlating the HRT to the TN removal efficiency of the systems, the explanatory power was even lower (r 2 < 0.07), except for Aakaer (r 2 = The differences between inlet and outlet for monthly average discharged-weighted nutrient concentrations were generally higher during the summer half year and highest for the hydrological year 2018/2019 (Figure 3). For Aakaer, TN removal efficiency fluctuated between 9% and 16% during 2016-2019, being lower compared to previous hydrological years ( Table 2). Removal efficiency of TP here was, on the other hand, the highest compared to the other SFW+WBs, and reached 79% in 2018/2019. For Ondrup, TN removal efficiency varied between 16% and 32%, while for TP it varied between −23% and 17%. For Ryaa, TN removal efficiency varied between 8% and 32%, with minimal effect in 2016/2017 and maximum in 2018/2019. Removal efficiency of TP varied between 0% and 44%, and was generally high, despite the low inflow and outflow TP concentrations.
The simple regressions carried out in this study demonstrated that HLR correlated inversely to TN removal efficiency in the SFW+WBs (Figure 4). However, in Aakaer (Figure 4a) a quadratic model explained nearly half of the variation in TN removal efficiency and indicated that an increase may be achieved after a certain HLR threshold. When correlating the HRT to the TN removal efficiency of the systems, the explanatory power was even lower (r 2 < 0.07), except for Aakaer (r 2 = 0.45 with slope = 2.34). In regards to TP, the variation of removal efficiency in the systems was weakly affected by both HLR (r 2 nearly zero) and HRT (r 2 < 0.12), except for Ryaa (r 2 = 0.33 with slope = 11.8 in relation to HLR; and r 2 = 0.50 with slope = −9.0 in relation to HRT), although these were not yet major explanatory variables.
Water 2020, 12, x FOR PEER REVIEW 8 of 15 0.45 with slope = 2.34). In regards to TP, the variation of removal efficiency in the systems was weakly affected by both HLR (r 2 nearly zero) and HRT (r 2 < 0.12), except for Ryaa (r 2 = 0.33 with slope = 11.8 in relation to HLR; and r 2 = 0.50 with slope = 9.0 in relation to HRT), although these were not yet major explanatory variables.
(a) (b) (c) The effect of TN load on TN removal rate was rather weak in all SFW+WBs (0.13 < r 2 < 0.33 with 0.06 < slope < 0.24) (Figure 5a-c). In contrast, TP load was generally a major explanatory variable for TP removal rate, especially in Ryaa (r 2 = 0.82 with slope = 0.81). In Ondrup, the latter correlation was weaker after disregarding two data points identified as outliers (r 2 = 0.33 with slope = 0.38) ( Figure  5e). Multiple regressions revealed that HLR and TN concentration explained little of the variation in TN removal rate and efficiency in the SFW+WB. Additionally, the individual effect of HLR was generally found significant in contrast to TN concentration (Table 3). A significant and relatively strong relationship was found in Aakaer in relation to TN removal rate and concentration. However, for TP the regressions demonstrated a strong relationship of both HLR and TP concentration in the variation of the TP removal rate, except in Ondrup. Correlations of the HLR with TP removal efficiency were weak, especially in Ondrup. Moreover, the individual effect of each variable was significant in the variation of TP removal rate for all systems, in which TP concentration was a particularly strong factor in Aakaer and Ondrup, as indicated by the relatively high coefficients. However, in relation to TP removal efficiency, significant correlations with individual variables were found in Aakaer and especially in Ryaa (i.e., the latter presented markedly high coefficients). The effect of TN load on TN removal rate was rather weak in all SFW+WBs (0.13 < r 2 < 0.33 with 0.06 < slope < 0.24) (Figure 5a-c). In contrast, TP load was generally a major explanatory variable for TP removal rate, especially in Ryaa (r 2 = 0.82 with slope = 0.81). In Ondrup, the latter correlation was weaker after disregarding two data points identified as outliers (r 2 = 0.33 with slope = 0.38) (Figure 5e).
Water 2020, 12, x FOR PEER REVIEW 8 of 15 0.45 with slope = 2.34). In regards to TP, the variation of removal efficiency in the systems was weakly affected by both HLR (r 2 nearly zero) and HRT (r 2 < 0.12), except for Ryaa (r 2 = 0.33 with slope = 11.8 in relation to HLR; and r 2 = 0.50 with slope = 9.0 in relation to HRT), although these were not yet major explanatory variables.
(a) (b) (c) The effect of TN load on TN removal rate was rather weak in all SFW+WBs (0.13 < r 2 < 0.33 with 0.06 < slope < 0.24) (Figure 5a-c). In contrast, TP load was generally a major explanatory variable for TP removal rate, especially in Ryaa (r 2 = 0.82 with slope = 0.81). In Ondrup, the latter correlation was weaker after disregarding two data points identified as outliers (r 2 = 0.33 with slope = 0.38) ( Figure  5e). Multiple regressions revealed that HLR and TN concentration explained little of the variation in TN removal rate and efficiency in the SFW+WB. Additionally, the individual effect of HLR was generally found significant in contrast to TN concentration (Table 3). A significant and relatively strong relationship was found in Aakaer in relation to TN removal rate and concentration. However, for TP the regressions demonstrated a strong relationship of both HLR and TP concentration in the variation of the TP removal rate, except in Ondrup. Correlations of the HLR with TP removal efficiency were weak, especially in Ondrup. Moreover, the individual effect of each variable was significant in the variation of TP removal rate for all systems, in which TP concentration was a particularly strong factor in Aakaer and Ondrup, as indicated by the relatively high coefficients. However, in relation to TP removal efficiency, significant correlations with individual variables were found in Aakaer and especially in Ryaa (i.e., the latter presented markedly high coefficients). Multiple regressions revealed that HLR and TN concentration explained little of the variation in TN removal rate and efficiency in the SFW+WB. Additionally, the individual effect of HLR was generally found significant in contrast to TN concentration (Table 3). A significant and relatively strong relationship was found in Aakaer in relation to TN removal rate and concentration. However, for TP the regressions demonstrated a strong relationship of both HLR and TP concentration in the variation of the TP removal rate, except in Ondrup. Correlations of the HLR with TP removal efficiency were weak, especially in Ondrup. Moreover, the individual effect of each variable was significant in the variation of TP removal rate for all systems, in which TP concentration was a particularly strong factor in Aakaer and Ondrup, as indicated by the relatively high coefficients. However, in relation to TP removal efficiency, significant correlations with individual variables were found in Aakaer and especially in Ryaa (i.e., the latter presented markedly high coefficients). Table 3. Multiple regression of the monthly (mo) effect of hydraulic loading rate (HLR) (m/mo) and total nitrogen or phosphorus concentration (TN and TP) (mg/L) as independent variables on the removal rate (g/m 2 /mo) and efficiency (%) of TN and TP for the SFW+WB in Aakaer, Ondrup, Ryaa, Gjol, and the SFW in Gjol.

SFW+WB vs. SFW in Gjol
The HLR variations of the SFW+WB and SFW in Gjol were similar throughout the study period (Figure 6a,b). However, the intensity of the HLR differed due to the Q ratio (~1:2.5) between the systems and their areas, corresponding to 1875 and 8100 m 2 in the SFW+WB and SFW, respectively. The highest HLRs were observed during the winter half year, while the lowest occurred during the summer half year. Peaks above 0.3 m/day occurred in the SFW+WB, while the highest HLR reached 0.2 m/day for the SFW. Differences in HRT between the two systems were more pronounced. Values of HRT for the SFW+WB were lower and equal to 3 and 7 days (average) for the winter and summer half year, respectively. Values of HRT for the SFW were higher and equal to 7 and 27 days (average) for the winter and summer half year, respectively, which are directly regulated by the volume of the systems, corresponding to 1263 and 6000 m 3 in the SFW+WB and WB, respectively.
Typical TN concentrations at the inlet of the two systems varied from 4 to 13 mg/L (Figure 6c,d). Outlet TN concentrations varied between 2 and 12 mg/L, and 2 and 13 mg/L for the SFW+WB and SFW, respectively. The differences between the inlet and outlet daily TN concentrations are higher during the summer half year. A total of 5 and 12 kg/ha corresponding to 27% and 26% of the TN load were removed during the monitored period by the SFW+WB and SFW, respectively. Common inlet TP concentrations varied between 0.37 and 2.21 mg/L. Outlet TP concentration varied between 0.19 and 3.51 mg/L, and 0.19 and 2.95 mg/L for the SFW+WB and SFW, respectively. Values of TP concentrations were higher in the summer half year and were more variable. Overall, there is not a substantial difference between inlet and outlet daily TP concentrations during the monitored period. However, a retention capacity of 0.05 mg/L corresponding to 6% of the TP load was measured in the SFW. The SFW+WB was capable of removing only small amounts of TP and an overall release of 0.11 mg/L corresponding to 24% of the TP load was measured.
The simple regressions demonstrated that the correlations were rather weak for both the SFW+WB (r 2 < 0.28) and SFW (r 2 < 0.17) located in Gjol. The strongest relationships were found between HLR and TN removal efficiency in the SFW+WB (r 2 = 0.33 with slope = 3.09 (x 2 ) and −34.15 (x)), and between HRT and TN removal efficiency in the SFW (r 2 = 0.29 with slope = 1.36). The multiple regressions demonstrated that a combined effect of HLR and TN concentration was not able to significantly explain the variations of TN removal rate and efficiency in both systems ( Table 3). The only exception was found for the SFW in relation to TN removal efficiency, which revealed not only that the combined effect largely explained the variation, but also that TN concentration was the only significant variable. For TP, on the other hand, HLR and TP concentration largely explained the variations of both TP removal rate and efficiency for all systems, especially in the SFW+WB, by which only the effect of TP concentration was significant. The simple regressions demonstrated that the parametric correlations were rather weak for both the SFW+WB (r 2 < 0.28) and SFW (r 2 < 0.17) located in Gjol. The strongest relationships were found between HLR and TN removal efficiency in the SFW+WB (r 2 = 0.33 with slope = 3.09 (x 2 ) and −34.15 (x)), and between HRT and TN removal efficiency in the SFW (r 2 = 0.29 with slope = 1.36). The multiple regressions demonstrated that a combined effect of HLR and TN concentration was not able to significantly explain the variations of TN removal rate and efficiency in both systems ( Table  3). The only exception was found for the SFW in relation to TN removal efficiency, which revealed not only that the combined effect largely explained the variation, but also that TN concentration was the only significant variable. For TP, on the other hand, HLR and TP concentration largely explained the variations of both TP removal rate and efficiency for all systems, especially in the SFW+WB, by which only the effect of TP concentration was significant.
The response of the SFW+WB and WB in terms of outlet nutrient concentrations is compared in Figure 7a. Despite the differences in the nutrient and hydraulic loads (Table 2), the removal response appeared almost identical for TN (r 2 = 0.96) and outlet concentrations followed the 1:1 line, indicating the condition of equal performance. Substantial differences and greater deviation from the 1:1 line were identified for the TP outlet concentrations (r 2 = 0.11), thus indicating a different performance for the two systems.
The contribution to the TN and TP removal/release (%) for two of the main compartments of the SFW+WB (SB and WB) is given in Figure 7b. The average daily TN removal is 33% and 5% for the WB and SB, respectively. Results show that that WB is more efficient compared to the SB in terms of TN removal, since most of the points are located above the 1:1 line indicating conditions of equal contribution to the total values presented in Table 2. The average daily TP retention is 23% and 6% for the WB and SB, respectively. Thus, the WB plays also an important role in removing TP. The response of the SFW+WB and SWF in terms of outlet nutrient concentrations is compared in Figure 7a. Despite the differences in the nutrient and hydraulic loads (Table 2), the removal response appeared almost identical for TN (r 2 = 0.96) and outlet concentrations followed the 1:1 line, indicating the condition of equal performance. Substantial differences and greater deviation from the 1:1 line were identified for the TP outlet concentrations (r 2 = 0.11), thus indicating a different performance for the two systems. The contribution to the TN and TP removal (%) for two of the main compartments of the SFW+WB (SB and WB) is given in Figure 7b. The average daily TN removal is 33% and 5% for the WB and SB, respectively. Results show that the WB is more efficient compared to the SB in terms of TN removal, since most of the points are located above the 1:1 line indicating conditions of equal contribution to the total values presented in Table 2. The average daily TP retention is 23% and 6% for the WB and SB, respectively. Thus, the WB also plays an important role in removing TP.

Nitrogen Removal Evaluation
Despite the low TN load, removal rates for the SFW+WBs estimated in this study fall within the removal rates found in previous studies in relation to woodchip-based systems [40,48]. A positive correlation was identified between HLR, TN removal rates, water flow, and seasonality which is also supported in the previous studies. Greenan et al. [33] and Robertson and Merkley [36] reported increasing N removal rates for increasing water flow. Hoffmann et al. [24] reported the highest N removal rates during the summer season for different WB designs. Removal efficiencies of TN in our investigated SFW+WBs varied between 8% and 32% during the years 2016-2019. Better performances were achieved during the previous hydrological years reaching the highest level of 51% at Ondrup. Addy et al. [48] carried out a meta-analysis on denitrifying bioreactors for nitrate removal and found that denitrifying bioreactors less than 13 months old had significantly higher nitrate removal rates than those being between 13 and 24, and >25 months old (p < 0.05). This may explain the high TN removal (%) for Aakjaer and Ondrup during the first monitored hydrological years (Table 2). Volumetric removal rates oscillate between 0.19 and 2.31 g/m 3 bioreactor/day. Schipper et al. [40] found that under field operating conditions, denitrifying systems could achieve a volumetric removal rate of 2 to 22 g/m 3 bioreactor/day with the lowest rate often associated with nitrate-limitations. Results from our single regression showed a negative correlation between the TN removal efficiency and HLR (Figure 4). The simple regression analysis further showed that TN load weakly explained the variation in removal rate (Figure 5a-c). This indicates the lack of nitrate-limiting conditions under these circumstances and the major role played by other driving factors on the process of N removal.
The meta-analysis of Addy et al. [48] showed that N removal for HRT less than 6 h was significantly lower than in beds with HRT from 6 to 20 h and more than 20 h (p < 0.05). Our results show an average HRT in the WB of less than 5 h during the winter half year, which is generally characterized by the highest load. Results of our simple and multiple regression analysis showed that hydrology (HLR and HRT) plays an important role in the TN removal ( Figure 4). However, hydrological parameters were not major predictors, indicating that other factors (e.g., temperature) may be dominant in regulating denitrification. Hoffmann and Kjaergaard [49] suggested incorporating hydraulic control components in the system design which are active during high flow events in order to increase the HRT. This is advisable to increase the nutrient removal during the winter half year and avoid unfavorable side effects such as methane or hydrogen sulphide emission from the system during the summer half year. Potential low oxygen concentrations in the WB outlet should also be addressed. This is particularly important for catchments with small vulnerable recipients, which are primarily fed by the WB outflow. In this situation, methods of aeration of the outflow should be applied.

Phosphorus Retention Evaluation
Recent studies proved TP retention in WBs, highlighting the ability of denitrifying-based systems to not only treat N. Gottschall et al. [43] estimated the performance of inline WBs for TP retention finding a median retention efficiency of 28%. Despite the low TP load, the average TP retention efficiency across the investigated SFW+WBs was 49% with peak retention performances >64%. The hydrological year 2018/2019 is not considered as it was particularly dry (Table 2). Volumetric retention rates oscillate between 0.01 and 0.15 g/m 3 bioreactor/day. A negative correlation was also identified between the low concentrations and the TP retention efficiency.
Choudhury et al. [44] demonstrated that almost all TP retention in woodchip filters (91%) was associated with PP of which 71% was removed by a sedimentation tank (12.3 m 3 ). This is in not in agreement with the results from Gjol (Figure 5b). Retention of TP in the SB appeared to be more occasional and higher (in %) while in the WB is more recurrent and lower. Hydrological parameters (HLR and HRT) have a large effect on P retention efficiency, especially in Ryaa. However, they are not major predictors for the TP retention. This differs with the results from a review on edge-of-field technologies in which the HRT was identified as a key factor regulating P retention [50]. Our results do not indicate saturation of the SB or clogging of the WB by P species throughout the monitored period. This indicates a lifespan of more than three years for the SFW+WB in terms of P retention capability under similar conditions.
The presence of a pumping system before the SFW+WBs may play a key role in defining the lifespan of the facility as it prevents a large amount of sediment and thus, PP to be treated, as observed in sedimentation basins preceding SFWs as pre-treatment systems [51]. This is the case for our paired systems, where drainage water is pumped from an open ditch into the SFW+WB and SFW. The absence of sediments and the controlled water supply at these facilities can guarantee longer and better functionality, thus avoiding overflow throughout the monitored period. In addition, temporal tile drainage water retention in pumped agricultural areas provides an instrument for the adaptation to climate changes. Targeted pump operations in fact allow keeping the water level at optimal conditions for plant growth and development, and for protection of sensitive areas. Results from a longer monitored period will allow a better comparison between the investigated SFW+WBs to define the pumping effect on the overall TP retention performance.

Conclusions
Agriculture is considered a primary source of diffuse pollution responsible for eutrophication followed by the loss of biodiversity in marine and freshwater ecosystems. Surface-flow constructed wetlands paired with woodchip bioreactors (SFW+WB) are promising nutrient edge-of-field technologies targeting nitrogen (N) and phosphorus (P) losses from agricultural drainage water. Results from this study showed that these systems provided good removal rates throughout the monitored period. In particular, removal of TN varied between 8% and 51% and lower removal rates were identified with increasing system age. Retention of TP ranged from 0% to 80% and was generally negative (nutrient release) during dry hydrological years. A negative correlation was identified between the nutrient concentrations and the removal rate. Additionally, the statistical analysis revealed that the hydraulic loading rate (HLR) and hydraulic retention time (HRT) were not major explanatory variables to the variation in nutrient removal rate.
A comparison of the performance of an SFW+WB and an SFW underlined the possibility of implementing the paired system in management plans at the national scale. Additionally, the WB is very effective in the removal of TN as it represents the carbon source activating the denitrification process. Overall, the proposed paired system showed a nutrient removal capacity comparable with that of SFW. The main advantage to implementing WBs is that smaller agricultural areas are needed for the system installation. Farther testing including analysis of concentration and proportion of nutrient particulate and dissolved forms, especially during highly variable and intense discharges, must be considered in future works. In addition, a cost-effective analysis is recommended for evaluating the synergetic effect of SFW+WB in comparison to standalone and currently used systems, their costs, and performance maximization.