Pesticide peak concentration reduction in a small vegetated treatment system controlled by chemograph shape

. Pesticides may impact aquatic organisms when entering water bodies. Measures for mitigation against pesticide 10 inputs include vegetated treatment systems (VTS). Some of these systems have very short hydraulic retention time (< 1 h) but nevertheless manage to effectively reduce peak concentrations of contaminants. We hypothesize that this is not solely the result of contaminant and VTS properties but also related to the shape of the contaminants input chemograph, i.e. its sensitivity to dispersion. In order to test this hypothesis we performed a cluster analysis on the chemographs of contaminants mobilized in response to rainfall events in a viticultural catchment and derived a measure for dispersion sensitivity which we included 15 into multiple linear regression analysis. The resulting measure was then incorporated into multiple linear regression models for description of contaminant mitigation in a VTS located at the catchment outlet. We found that the mobilization clusters reflected the source areas of the contaminants and that dispersion sensitivity was the dominant explanatory variable for the reduction of contaminant peak concentration. These findings imply that mitigation of the toxicological risk in VTS may be stronger for compounds with pronounced peaks than for those abundantly available in a catchment.


Introduction
Use of pesticide is beneficial for agricultural productivity. When pesticides reach surface water bodies, they threaten aquatic organisms (Zubrod et al., 2019). Effects of pesticides on aquatic ecosystems include a reduction of species richness of invertebrates (Beketov et al., 2013) as well as microorganisms (Fernández et al., 2015). Unintended export of pesticides from the application site to water bodies can happen in particulate form via erosion (Oliver et al., 2012;Taghavi et al., 2011) or in 25 dissolved form via surface runoff, drainage pipes, spray drift or leakage to groundwater (Reichenberger et al., 2007) and subsequent exfiltration.
In the environment pesticides are subject to degradation by both abiotic (e.g. hydrolysis, photolysis) and biotic (e.g. plant metabolism, microbial degradation) processes (Fenner et al., 2013). If degradation is incomplete, pesticides form transformation products (TPs) which may be equally or more mobile, persistent and toxic than their PCs (Hensen et al., 2020). 30 from that of their parent compounds (PCs) in terms of both source areas and export pathways as the formation of TPs may happen on larger time scales and TPs mostly have different physicochemical properties than their PCs .
The specific transformation and further degradation of a contaminant largely depends on the interplay of the contaminant's 35 mobility and degradability as well as site characteristics (Gassmann et al., 2015). Both mobility and degradability can vary over multiple orders of magnitude for different contaminants. However, water and soil half-lives are at least in the order of several days or weeks for most pesticides (Lewis et al., 2016).
Pesticide pollution can be mitigated by vegetated treatment systems (VTS) located between source areas and receiving water bodies (Vymazal and Březinová, 2015;Stehle et al., 2011). Such systems temporally retain polluted waters and thus provide 40 space, time and favorable conditions for degradation processes to happen. VTSs studied in literature include very different types of systems (Lange et al., 2011), including vegetated ditches or detention ponds with hydraulic retention times (HRT) ranging in the order of minutes to several hours (Bundschuh et al., 2016;Elsaesser et al., 2011) or constructed wetlands in which HRT may reach several weeks, when operated in batch mode (Tournebize et al., 2017;Maillard et al., 2016). Pesticide mitigation has also been observed in other systems in which water is temporally stored such as barrage fish ponds (Gaillard et 45 al., 2016) or rice fields (Moore et al., 2018). While contaminant mass loss (RM) is mainly observed in systems with longer HRT, systems with short HRT have been shown to efficiently reduce peak concentration (RC) and the associated acute toxicity (Bundschuh et al., 2016;Elsaesser et al., 2011;Stehle et al., 2011).
During longitudinal transport in streams, peak concentration reduction does not necessarily involve degradation but can be the result of enhanced dispersion due to the presence of obstacles such as plants (Elsaesser et al., 2011) and removal from the 50 water phase by sorption. Mitigation properties therefore constantly change due to wetland succession . However, following the concept of advective-dispersive transport (Fischer et al., 1979), the effect of dispersion on a concentration signal in terms of peak concentration reduction does not only depend on the dispersion imposed on that signal but also on the shape of the signal. Peak concentration reduction will be stronger for a signal with pronounced peak and low background than for a signal with small peak and high background if both signals are exposed to the same dispersion. The 55 influence of this effect on variability in contaminant mitigation and its implications for VTS efficiencies have not been systematically investigated so far.
In this study, we assess the mobilization of six organic contaminants in response to rainfall in a viticultural catchment and the influence of the contaminant concentration signal on mitigation efficiency in a VTS located at the catchments outlet. We hypothesize that mitigation efficiency does not only depend on properties of the contaminant (e.g. sorption affinity, water and 60 soil half-lives) or the treatment system (e.g. retention time, plant coverage), but also on the shape of the concentration signal mobilized in the catchment. To test this hypothesis, we search for patterns in contaminant flush signals by performing a cluster analysis on flow-triggered monitoring data, integrate the results into two multiple linear regression models for peak concentration reduction and mass removal, respectively, and evaluate the relative importance of the model variables.

2
Material and methods 65

Study site
The study site ( Fig. 1) is located inside a flood detention basin in a 1.8 km² headwater catchment, southwest Germany.
Catchment elevation ranges from 213 to 378 m.a.s.l.. Mean precipitation was 800 mm a -1 between 2009 and 2018, mean air temperature was 11.3 °C. Soils mainly consist of calcaric regosol which formed on Aeolian loess and have a typical grain size distribution of 10 % sand, 80 % silt and 10 % clay. Most of the area is dedicated to viticulture on large artificial terraces (71 %), 70 while croplands occupy the main valley (20 %). Forest only accounts for a small portion (9 %) and is limited to the most elevated part of the catchment. The elevated vineyard terraces are subject to frequent fungicide application, while herbicides are applied to the cropland in the flat valleys. Large parts of the catchment are drained by a dense sub-surface pipe network with a total length of about 9 km which causes fast downstream transport of storm water as well as dissolved and suspended material (Gassmann et al., 2012). 75 Inside the detention basin, a 258 m² vegetated surface flow constructed wetland and a 105 m² retention pond (maximum depth 1.5 m) are connected in series parallel to the course of the Löchernbach stream. During baseflow conditions a small dam diverts all flow through the vegetated treatment systems. In the case of large discharge events, the treatment systems are bypassed.
The wetland is in operation since 2010 and its succession was studied by Schuetz et al. (2012). The pond was added to the system in January 2016. Both segments are sealed by an impermeable clay layer that prevents infiltration. From previous 80 studies, it is known that HRT in the system ranges from about one hour during flood events to more than a days during extreme low flow.

Monitoring setup
Stream flow was measured every minute between April 2016 and September 2017 at two gauges about 200 m upstream of the treatment system (G1) and at its outlet (G2). Water levels at G1 were recorded inside a 1.37 m standard H-flume (Bos, 1989) by means of a pressure transducer (Decagon CTD10) and related to discharge using a standard rating curve. At G2 water levels were measured in a rectangular cross-section 2 meters ahead of the detention basin outlet by a radar gauge (Vegapuls 61). The 90 corresponding rating curve for G2 accounted for complete submergence of the control gate valve (Peter, 2005).
A total of 15 sampling campaigns were performed (Fig. 2), ten of which assessed rainfall induced contaminant mobilization and five represented stationary flow. During stationary conditions, grab samples were taken at G1 and G2. Duplicates were produced by splitting the sample into 2 brown glass bottles with a volume of 1 l each. During transient flow conditions pesticide sampling was automated. When the upper gauge registered a water level increase of more than 3 cm/h, an automatic sampler 95 (ISCO 3700) started to fill pairs of 900 ml glass bottles at 0, 0.5, 1, 2, 6, and 12 hours after activation. A second automatic sampler was launched at G2 following the same sampling scheme with a time lag of one hour to account for transit between G1 and G2. All samples were recovered from the study site within 24 h after sampling and cooled until analysis.

Analytical methods and pesticide analysis
The target compounds included the two fungicides boscalid and penconazole, the two herbicides metazachlor and flufenacet, and the two TPs metazachlor sulfonic acid (met-ESA) and metazachlor oxalic acid (met-OA). Selected physicochemical properties of the target compounds are listed in Table 1. According to the Pesticide Properties Data Base (Lewis et al., 2016) https://doi.org/10.5194/hess-2020-228 Preprint. Discussion started: 2 July 2020 c Author(s) 2020. CC BY 4.0 License.
the contaminants can be classified as low (boscalid) to moderately soluble in water. Mobility ranges from very mobile (TPs) 105 to slightly mobile (fungicides). The fungicides are considered moderately fast degradable in the water phase and persistent in soils, while the herbicides are considered stable in the water phase and non-persistent in soils. TPs of metazachlor are considerably more persistent in soil than their PC. The fungicides are considered stable with respect to hydrolysis but degradable via photolysis, while the herbicides are stable regarding both.

Identification of patterns in input concentrations
Identification of patterns in input concentration was done by k-means cluster analysis of the contaminant input concentration sequences recorded at G1. This method was first applied in signal processing and became a popular tool in data mining where it is used for pattern recognition and also found its way into hydrology, e.g. k-means clustering has recently been applied to 135 stream nitrate time series data (Aubert and Breuer, 2016). The k-means approach partitions the elements of a dataset into a predefined number k of clusters by attributing the elements to the cluster with the nearest cluster center. Optimal clustering is achieved by iteratively updating cluster centers and minimizing within cluster variance.
A total of 58 concentration sequences was included in the analysis, consisting of 10 sequences per target compound, except for flufenacet which did not exceed LOQ in two events. Prior to cluster analysis, data was normalized by the maximum of 140 each concentration sequence to guarantee that clustering was done by signal shape, not by total concentration. The analysis was done in the software R (R Core Team, 2019) using an algorithm by Hartigan and Wong (1979). We tested clustering for k ranging between 2 and 10, the final number was determined by both visual inspection of the clusters and assessment of explanatory benefit per additional cluster (elbow method). We found that k=4 resulted in the best partition.
In order to integrate the results of the clustering into further analysis, we developed the dispersion sensitivity index iDS. 145 Following the concept of advective-dispersive solute transport, relative peak concentration reduction can be expected to be higher for concentration signals with clearly defined peaks than for those with relatively flat peaks compared to the background.
The dispersion sensitivity index was correspondingly defined as the relative portion of peak concentration (Ĉ in ) susceptible to dispersion according to Eq. 1: where , is the concentration in the last sample collected, considered as a reference for post-peak concentration. 150
where ̂ and ̂ are peak concentration registered at the inlet and outlet sampling points. RM was calculated analogously 155 from the input (Min) and output contaminant mass (Mout) as shown in Eq. 3: Contaminant masses were calculated from discharge at G1 and G2 and linearly interpolated contaminant concentrations. For comparability with automated sampling, contaminant masses during stationary flow were referred to a period of 12 h during which constant concentrations and flow were assumed. As water level data from G2 showed evidence for inaccuracy during low flows due to the constant cross-section, we assumed that the VTS was in equilibrium during stationary flow conditions 160 and used flow from G1 for calculation of both Min and Mout. Following this procedure, mass removal assessment is limited to the liquid phase as contaminants adsorbed to sediments are not taken into account.

Identification of influential variables
In order to identify influential variables for contaminant mitigation we constructed two separate multiple linear regression models for RC and RM and attributed relative importances to the model variables. As TP peak concentration data turned out to 165 be affected by low sampling frequency during late stages of the sampling procedure, we decided to focus on PCs and exclude TPs. This decision reduced the number of data points but prevented the model from being affected by the high variability in TP mitigation efficiencies resulting from imperfect sample coverage. Model construction was done as follows. First, an initial set of potential model variables was selected based on literature and practical considerations for each model. Then all-subsets regression (Lumley, 2017) was performed on the initial set of variables, i.e. all possible numbers of variables and variable 170 combinations were tested and evaluated in terms of explanatory power as expressed by adjusted R².
The initial variable sets of the two models are represented by the light columns in the background of Fig. 6. They included physicochemical contaminant properties such as the organic carbon adsorption coefficient (KOC) and water solubility (Solub).
In contrast to other studies (Bundschuh et al., 2016;Stehle et al., 2011), we decided to not include the degradation half-lives listed in Table 1 as even the lowest of these values (Water T50 of penconazole is reported to be 2 days) would not have produced 175 measurable effects in the short retention time of approx. one hour during event conditions. Hydraulic conditions were represented by mean discharge ( ̅ ) and mean hydraulic retention time (HRT). Variables considered relevant for only one of the models were input peak concentration (̂) and dispersion sensitivity (iDS) in the RC model as well as input mass ( ) and the water balance error eW which was included in the RM model. A variable identified as important by other studies (Stehle et al., 2011;Bundschuh et al., 2016;Elsaesser et al., 2011) but not included in our study is vegetation coverage. VTS properties 180 such as vegetation coverage may be essential for comparison of different systems (Lange et al., 2011;Stang et al., 2014) or states of systems . However, this was not the main purpose of our study.The data set was checked for multi-collinearity of independent variables by calculating Pearson correlation coefficients and variance inflation factors (VIF).
Strongest correlation was found between ̅ and HRT (R=0.85). VIF for these variables were between 3 and 5. This is below the standard threshold of 10, but as Hair (2010) points out that in models with limited data, comparatively small VIF values 185 https://doi.org/10.5194/hess-2020-228 Preprint. Discussion started: 2 July 2020 c Author(s) 2020. CC BY 4.0 License. may indicate substantial collinearity, we decided to keep the variables that was stronger correlated with the mitigation rates ( ̅ ) and exclude the other one (HRT) from the model building routine. The exclusion of HRT, however, only marginally changed the model output.
Data points were considered outliers and removed from the data set, if at least one of the following condition was met: (1) Cook's distance exceeded one, indicating particularly influential points, or (2) RC or RM, in the respective models, deviated 190 from the mean by more than 2 standard deviations. The number of observations in the data sets passed to all-subset regression was 36 for the RC model and 37 for the RM model. Relative importance of the predictors in the resulting model was determined by decomposing total R² into non-negative contributions based on the corresponding sum of squares (Grömping, 2006). As attribution of explanatory power to predictors can depend on the order in which they are added to the model, we applied the LMG measure (Lindemann et al., 1980) which averaged predictor importance over all possible orderings. 195

3
Results and discussion

Contaminant mobilization patterns
Clustering of input contaminant breakthrough curves by shape similarity resulted in an optimal number of four clusters (Fig.   3). These clusters differed in terms of shape, peak timing, dispersion sensitivity and contaminant composition. Cluster A mainly represented the fungicides (boscalid and penconazole) was characterized by a short median peak arrival time (Tpeak) of 200 1 h and an elevated shoulder, i.e. concentration often had not dropped back to zero by the end of the sampling procedure. Due to the pronounced tailing, the dispersion sensitivity index reached intermediate values. Cluster B was of mixed composition and showed a similarly quick increase (Tpeak=1 h). Cluster B differed from cluster A by a lower tailing and corresponding higher dispersion sensitivity. Cluster C was dominated by the TPs of metazachlor and was overall less concise but clearly distinguishable from the other clusters by the late peak (Tpeak=6 h). In cluster C peak and late-time concentration were often 205 similar which resulted in very low dispersion sensitivity. Cluster D consisted mainly of the herbicides (metazachlor and flufenacet) and was similar to cluster B but peak concentration was reached about one hour later (Tpeak=2 h). With the exception of cluster B which rather represented similar events (event 1 and event 4 in Fig. 2), overall clustering was controlled by similar behavior of contaminant groups.
We interpret these findings as an indication of different source areas and input transport pathways. In the study catchment 210 fungicides (cluster A) are in fact mainly applied to the elevated vineyard terraces and likely to quickly reach the stream via surface runoff and the drainage network. Herbicides (cluster D), in contrast, are rather applied to croplands in the flat valley areas and transported to the stream more slowly due to lower terrain slope. Rising water levels in the stream may further reduce the hydraulic gradient in drainage pipes and thereby reduce flow velocity. The TPs (cluster C) have the same geographical source area as their PC but are formed in the soil and seem to be transported comparatively slowly to the drainage network via 215 matrix leaching. This interpretation is in line with findings from a field study by Doppler et al. (2012) who showed that a combination of macropores and tile drains may act as a quick export pathway of pesticides to the stream. Slower export pathways of TPs compared to PCs were also reported by Gassmann et al. (2013) who found that soil matrix leaching rates to https://doi.org/10.5194/hess-2020-228 Preprint. Discussion started: 2 July 2020 c Author(s) 2020. CC BY 4.0 License. tile drains played a bigger role for TPs than for PCs in a Swiss catchment. This finding was attributed to higher mobility and higher amounts of initial residues of TPs compared to PCs in the soil matrix. These explications also seem plausible in the 220 catchment investigated in this study.
One ambiguous aspect to k-means cluster analysis is the selection of the number of clusters. While we found k=4 to result in interpretable clusters, other numbers of clusters would be equally justifiable based on explanatory power gained per additional cluster. In order to make this somewhat arbitrary decision more transparent, we also checked changes in cluster composition if three or five clusters were allowed. We found that in both cases the fungicides and TPs were still largely separated, while 225 the remaining clusters represented intermediate properties, i.e. the overall interpretation of the clusters would not change. In the present study the function of the cluster analysis was mainly to illustrate differences in contaminant behavior, however, our results indicate a potential for its application for data exploration in catchments where less prior knowledge on catchment hydrology is available.

Input concentrations and loads
Input concentrations (G1) differed clearly under the different flow conditions. During stationary flow, boscalid and the TPs of 235 metazachlor were detected in almost all samples, while penconazole, metazachlor and flufenacet were only found occasionally.
Concentrations were usually in the order of tens of nanograms per litre (Fig. 4). During transient flow, all contaminants were found in all samples. Concentrations were generally much higher than during stationary flow and the concentration difference https://doi.org/10.5194/hess-2020-228 Preprint. Discussion started: 2 July 2020 c Author(s) 2020. CC BY 4.0 License. depended on the contaminant. Compared to stationary flow, median concentrations of boscalid increased by a factor of 48, while concentrations of met-ESA and met-OA only increased by a factor of 3 and 5, respectively. Medians of penconazole, 240 metazachlor and flufenacet were not calculated due to lack of data during stationary flow. The increase in concentrations translated into an even stronger increase in input load. During stationary flow, input mass related to a 12 hour period usually was in the order of tens of micrograms. During transient flow, input mass increased by about 4 orders of magnitude in the case of boscalid and about 2 orders of magnitude in the cases of met-ESA and met-OA. In extreme cases, several grams of boscalid (and metazachlor in one exceptional case) entered the treatment system during a single event. 245 Concentrations detectable at the catchments may depend on many factors including stream size (Lorenz et al., 2017) and landuse which complicates direct comparison. Peak concentrations of fungicides detected in this study, however, were in the range of those detected by Bundschuh et al. (2016) in different streams draining viticultural catchments in Southwest Germany. The fact that highest peak concentrations were detected for the fungicides reflects both the land-use distribution in the catchment (higher percentage of vineyards compared to cropland) and the different application practices (higher application frequencies 250 of fungicides). The increase in PC to TP concentration ratio from stationary to transient flow conditions indicates that TPs are exported via a permanently active pathway such as groundwater. In contrast, PCs are rather transported via periodically active pathways that respond quickly to rainfall. This interpretation is in line with the patterns identified in section 3.1 and findings of Gassmann et al. (2013).  The finding that RM was similar to RC is surprising at first sight because degradation of contaminants usually happens on larger 270 time scales in relation to HRT in the studied VTS or comparable systems (Elsaesser et al., 2011). Assessment of mass loss in our study, however, was limited to the liquid phase. Contaminant mass loss could therefore be the result of e.g. incorporation of contaminants by plants or adsorption to organic matter inside the VTS (Stang et al., 2014) or temporal trapping of highly contaminated portions of flow. Temporal trapping may have occurred as water flushed into remote areas of the detention basin where it was retained for durations exceeding the end of our event sampling scheme. Small, disconnected depressions filled 275 with water were in fact observed in the field after flood events. If such areas are temporally flooded, degradation conditions may be similar to those in VTS operated in batch mode, i.e. alternating oxic-anoxic conditions enhanced degradation as found by Maillard et al. (2016). Considering the different mobilization behavior of the contaminants, mass removal due to trapping may be different for the contaminants depending on which portion of flow was intercepted.
Experimental uncertainties may be associated to different discharge measurement methods at the two gauges and resulting 280 underestimation of contaminant mass. This in turn, may have caused overestimation of mass loss. However, relative water balance errors were usually smaller than RM. In order to check for a systematic relationship of water balance and RM, we included water balance as a potentially influential variable into the regression model building procedure. The sampling scheme seemed not to be the major source of error for PC mass, as PC concentrations had usually recessed to background levels when sampling ended (Fig. 3). This, again, was different for the TPs that arrived later and whose mass balances were therefore not 285 equally well captured by the sampling scheme.

Dispersion sensitivity dominates peak concentration reduction
The linear regression models explained 66 and 51 percent of variance in RC and RM, respectively (Fig. 6). This means that the models generally described concentration reduction more accurately than mass removal. The models also differed in terms of importance attributed to the variables. In terms of relative variable importance, the RC model of parent compounds was 295 dominated by the effects of iDS which accounted for about half of the explained variance. A secondary significant variable was discharge, contributions from input concentration and KOC were not significant. Mass removal was best described by a combination of discharge and input mass. The latter, however, was not significant, nor was KOC. The water balance error was excluded by the all subset regression model building procedure because it only marginally improved explanatory power of the RM model. 300 According to the RC model, RC was positively related to iDS, and negatively to mean discharge and input concentration. These results may be considered physically reasonable as high peak concentrations are likely to coincide with pronounced tailings and therefore low dispersion sensitivity. The importance of discharge in our case is considered equivalent to the importance of HRT found in other studies (Bundschuh et al., 2016;Stehle et al., 2011) as these variables were subject to strong negative correlation and are physically related. We therefore regard these results as a confirmation of our hypothesis that the input 305 signal shape should be relevant for RC. https://doi.org/10.5194/hess-2020-228 Preprint. Discussion started: 2 July 2020 c Author(s) 2020. CC BY 4.0 License.
In contrast to other studies, we did not find a significant relationship of RC with KOC (Stehle et al., 2011) or solubility (Bundschuh et al., 2016). In fact, the small contribution to explanatory power of KOC is based on a negative relationship, which contradicts the general opinion in literature (Vymazal and Březinová, 2015) and is probably an artifact resulting from the low number (n=4) of compounds and thus different KOC values included in the model. While iDS clearly improved the RC model, 310 no such specific variable was found for the RM model. The lower explanatory power of the RM model indicates that some important explanatory variable may be missing, possibly because temporal trapping of flow portions was not represented by the selected model variables. Conclusions In this study, we have shown that chemograph shapes of mobilized contaminants can be attributed to contaminant groups and 320 their source areas. Both comparison of absolute reduction rates and regression analysis suggest that the shape of the input signal may play an important role for peak concentration reduction in VTSs with short HRT. We therefore recommend that this factor should be considered for future assessment of VTS functionality and ecotoxicicty based risk assessment. As dispersion sensitivity is inherent in the concept of advective-dispersive transport, we consider these findings to be transferable to other solutes mobilized in catchments during storm events. Our findings imply that peak concentration reduction in VTSs 325 may generally be more efficient for compounds whose chemographs are characterized by sharp peaks and low backgrounds rather than for those ubiquitous in a catchment.

Data availability
Contaminant data used in this study is available as supplementary material.