Contribution of preferential flow to tile drainage varies spatially and temporally

Tile drainage of agricultural fields is a conduit for nutrient losses. Preferential flow in the soil can more directly connect surface runoff with tile drainage compared with matrix flow, which may also increase P losses. In this study, water temperature was monitored in surface runoff and tile drainage and electrical conductivity (EC) was measured in tile drainage at two sites in southern Ontario with different soil types (i.e., clay and loam). These data were used to estimate the percentage of preferential flow in tile drainage based on end member mixing. Estimates using temperature were compared with estimates using EC, and both were evaluated across seasons and hydrographs and compared with P concentration and load data. There was strong correlation (r = .83) between estimates of preferential flow using the two methods, but due to variability in surface temperatures, EC provided a less flashy estimate for preferential flow (Durbin–Watson statistics of 0.34 for temperature and 0.09 for EC). Preferential flow accounted for a higher percentage of tile drainage flow in clay soil than loam, but percentages were not significantly different between seasons or timing within events. Phosphorus concentrations and loads were weakly correlated with preferential flow, suggesting that P transport was influenced by other factors as well. Although further work is necessary to calibrate these methods for estimating preferential flow from continuously monitored temperature and EC, this technique can be applied to already collected data to model and test posited explanations of observed phenomena in P, other nutrients, and water transport from tile‐drained agricultural land.

Numerous in situ studies have linked preferential or macropore flow with P transport in tile drained landscapes (Heathwaite & Dils, 2000;Simard, Beauchemin, & Haygarth, 2000;Stamm et al., 1998;Williams et al., 2016). This is due in part to the decreased contact time and contact surface area with the soil, which decreases opportunities for the sorption of soluble reactive P (SRP), as well as saturation of sorption sites along macropores (Heathwaite & Dils, 2000;Stamm et al., 1998). High velocities in macropores also transport higher sediment loads than flow through the soil matrix because of a decrease in tortuosity (Heathwaite & Dils, 2000). These suspended sediments carry higher total P (TP) loads bound to and within aggregates compared with the dissolved fraction (Chapman, Foster, Lees, Hodgkinson, & Jackson, 2001). In a comparison of surface and tile concentrations of SRP and TP, many studies have found greater concentrations of both forms of P in surface flow Stamm et al., 1998). Since macropores provide a more direct connection between drainage tiles and the surface, preferential flow likely carries higher concentrations of P more characteristic of surface flow

Core Ideas
• Preferential flow in tile drainage was estimated using EC and temperatures. • Temperature and EC produced similar estimates of preferential flow in clay and loam. • The EC was a more consistent estimator of preferential flow than temperature. • Preferential flow was not significantly different between seasons or stages of events. • Phosphorus load correlated with preferential flow but was influenced by other factors.
compared with matrix flow through the soil (Dominguez, Bohlen, & Parmelee, 2004;Katterer, Schmied, Abbaspour, & Schulin, 2001;Singh & Kanwar, 1991). As tile drains contain a mix of this directly connected, high-P water and the lower-P matrix flow, average P concentrations in tile drains would be lower than in surface flow. However, quantifying preferential flow at large scales, both spatial and temporal, has proven difficult (Allaire, Roulier, & Cessna, 2009). Most studies to date have used tracers applied to small plots before an event (Grant, Macrae, Rezanezhad., & Lam, 2019;Stamm et al., 2002). Although this is effective at measuring residence time for an individual event, it is difficult to measure successive events with this method because applied tracers persist in the soil and there are a finite number of conservative tracers that can be used to measure subsequent events . This method also requires application of the tracers directly before events, which has logistical and budgetary constraints and changes antecedent conditions (Grant, Macrae, Rezanezhad & Lam, 2019). Measuring preferential flow using naturally occurring tracers still requires sample collection and therefore limits the temporal resolution of the data (Williams et al., 2016). Consequently, these factors have limited efforts to understand preferential flow contribution at annual and regional scales.
Measuring preferential flow with continuous sensor data has the potential to address the issues resulting from tracer applications, timing complications, and successive events. Sensors also allow for higher frequency data collection than discrete sample collection (Smith & Capel, 2018). Using sensors to measure preferential flow requires identifying naturally occurring, regularly replenished, relatively inert, conservative, and detectable tracers that have distinct end members. Temperature has been used to estimate water mixing in streams as described by Arrigoni et al. (2008). Indeed, studies have shown that temperature can be an effective passive tracer for stormwater in karst (O'Driscoll & DeWalle, 2006;Zajicek et al., 2011) and urban (Somers, Bernhardt, McGlynn, & Urban, 2016) watersheds, as well as agricultural tile drainage (Pluer, Morris, Walter, & Geohring, 2019). There is more extensive evidence of using electrical conductivity (EC) for hydrograph separation of rainfall input from soil and groundwater contributions in urban (Pellerin, Wollheim, Feng, & Vorosmarty, 2008), alpine (Cano-Paoli, Chiogna, & Bellin, 2019;Laudon & Slaymaker, 1997), agricultural (Michaud, Poirier, & Whalen, 2018), and mixed-use (Okello et al., 2018) watersheds. Indeed, investigations of flow paths using EC in field-scale tile-drained agriculture have successfully differentiated preferential flow from slow flow paths (Chikhaoui, Madramootoo, Eastman, & Michaud, 2008;Smith & Capel, 2018). In this study, we evaluated the use of temperature and EC as tracers for measuring preferential flow at the field scale. Here, we consider preferential flow to be water moving through the soil profile with very minimal interaction with the soil itself (Beven, 2018;Flühler, Durner, & Flury, 1996). Although many studies, including some of those mentioned above, delineate rainfall contribution, new flow, or event flow, we differentiated between event flow moving through macropores from event flow that moves through the soil matrix. We also compared the estimates of preferential flow using both techniques using high-frequency data collected during baseflow conditions and tile flow events over 4 yr through all seasons, and at two sites with different soil types. We then compared estimates of preferential flow volume with measurements of SRP and TP loads in tile drainage to understand the relationships between flow paths and P transport.

Study sites
Two sites in the Lake Erie watershed that have been used for ongoing agricultural runoff research Plach et al., 2019;Van Esbroeck, Macrae, Brunke, & McKague, 2016) were selected for further monitoring of preferential flow (Figure 1 of 900 mm (10% as snow) and 1,024 mm (30% as snow) (ECCC, 2019). Both sites were tile drained, with contributing areas of 7.4 and 8.7 ha for the clay and loam sites, respectively. The clay site had slightly shallower drains and narrower spacing (0.7 and 10.7 m, respectively) compared with the loam site (0.9 and 13.9 m, respectively). The sites had similar soil Olsen P values of 13 and 12 mg kg −1 at the clay and loam site, respectively. Both fields had similar crop rotations (corn [Zea mays L.]-soy [Glycine max (L.) Merr.]wheat [Triticum aestivum L.]) and fertilization (inorganic fertilizer with P three times during the study period, with application rates averaging 53 kg P ha −1 at the clay site and 161 kg P ha −1 at the loam site). Both sites also experienced shallow conservation tillage during the study period (clay every fall, loam in one fall), which likely disrupted macropores, but only for a short period (Williams et al., 2016). Further details on climate and management practices at the two sites are provided by Plach et al. (2019).

Data collection
Meteorological data were collected at each site with an automated weather station (Onset Corporation HOBO Weather Station). Stations were equipped with a HOBO U30 GSM data logger (logging at 15-min intervals) connected to a tipping bucket rain gauge (0.2-mm SRGB-M002 rainfall smart sensor).
Tile flow data and temperatures of surface runoff and tile drainage were collected for 26 storm and melt events at the clay site and 20 events at the loam site, all of which included both tile and surface flow ( Figure 1). Tile discharge was measured using simultaneous velocity and depth measurements with Hach Flo-Tote 3 sensors and FL900 data loggers (Hach, flow accuracy of ±5%) at 15-min intervals. HOBO level loggers (Onset, depth accuracy of ±0.1%) were used to verify flow depth in the tile drains, especially at low flow (depth < 2 cm), which was converted to flow based on stage-discharge relationships developed from periodic flow measurements of tile outlet during baseflow. Both the flow and depth sensors were secured to a mounting ring expanded inside the tile drains. The HOBO loggers were also used to measure water temperature continuously in the tile water, which covered all events. At each site, surface runoff was routed into a culvert or outlet pipe at the edge of the field to convey runoff for that single field separately from the tile drainage. A second HOBO logger was secured to a similar mounting ring within the culvert, which recorded the depth and temperature of surface runoff from each site for each event. ES-2 EC sensors (Decagon, EC accuracy of ±10%) were placed in the tile drains with the other sensors to measure EC of tile flow. This recorded data at 15-min intervals. Conductivity was available for 15 events at the clay site and 13 events at the loam site, all of which had corresponding temperature measurements.
For a subset of events, water samples were collected for water chemistry analysis. Only water chemistry from events with at least six samples during the hydrograph were used for statistical analyses. Automated water samplers (Teledyne ISCO 6712) were installed to collect from the tile drains and surface runoff collection pipes and were activated by water depth. Time-paced sampling intervals were used, adjusted seasonally, to capture entire event hydrographs (Van Esbroeck et al., 2016). Samples were removed from the ISCOs within 24 h and divided into aliquots that were either acidified (0.2% final H 2 SO 4 concentration) or filtered with 0.45-µm cellulose acetate filters and frozen prior to analysis. Filtered samples were analyzed for major anions (Cl − , CO 3 2− , NO 3 − , and SO 4 2− ) and cations (Ca + , K + , Mg 2+ , Na + , and NH 4 + ) using a Dionex ICS3000 (Thermo Fisher Scientific). Ions were analyzed for 72 samples during 12 events with both temperature and EC measurements at the clay site and 42 samples during three events at the loam site. Acidified aliquots were digested using the Kjeldahl method (Bowman & Delfino, 1982). Digested and filtered aliquots were analyzed for TP and SRP, respectively, using colorimetric analysis (Seal Analytical Methods no. G-188-97 and G-175-96 Rev. 13). Phosphorus was analyzed for 21 and 9 events (158 and 68 samples) with temperature data at the clay and loam sites, respectively. This included 14 and 5 events with P data (92 and 42 samples) with both temperature and EC data at the clay and loam sites, respectively. Concentrations measured below the detection limit were assumed to be half of the detection limit for that particular parameter and analytical method for statistical analyses.

Statistical analysis
The R software environment version 3.5.1 was used for calculations and statistical analysis. Flow, temperature, and EC data were averaged at an hourly interval to better fit with chemistry data. Events were delineated based on tile drain flow using the EcoHydRology package (Fuka, Walter, Archibald, Steenhuis, & Easton, 2018). This method separates baseflow from event flow by automatically creating recession curves as described by Nathan and McMahon (1990). The magnitude and timing of maximum flow was identified for each event. The period between the start of the event and the peak flow was considered the rising limb and the period between the peak flow and the event end was considered to be the falling limb. The periods 12 h before and 24 h after events were included in the calculation of baseflow, but not as part of the event flow or P load. Estimated matrix flow temperature was calculated as a linear interpolation between tile temperatures at the start and end of each event. Preferential flow temperature was assumed to be directly connected with surface water and therefore equal to the temperature of surface flow while there was surface flow. We assumed that temperature remained constant after surface flow stopped. Although variable value is a violation of the original assumption of fixed end members, it has been shown effective by Burns et al. (2001). Preferential flow EC was assumed to be 100 µS m −1 based on an average of periodic grab samples of surface water during events (n = 12, EC = 99.9 ± 4.8 µS m −1 ). Matrix flow EC was calculated using the same method as matrix flow temperature. These assumptions reflected a common approach for estimating tracer concentrations in preferential flow (Steenhuis, Boll, Shalit, Selker, & Merwin, 1994;Stone & Wilson, 2006). For both the temperature and EC methods, the percentage of flow at each timestep that was considered preferential flow was estimated as its value between the two end members (i.e., matrix and preferential). Thus, the quantity of preferential flow was calculated using the equation or EC [µS m −1 ]) of the matrix end member, E p is the value of the surface flow, which is assumed to describe the preferential flow, and x t is the observed value in the tile during the event (Stone & Wilson, 2006). Percentages <0% and >100% were assumed to be 0 and 100%, respectively. An assumption of end member mixing analysis is that end members show distinct signatures (Bansah & Ali, 2017;Hooper, Christophersen, & Peters, 1990). This assumption was tested between baseflow tile drainage and surface runoff for temperature, EC, and major ions. For each event, Mann-Whitney tests were used to compare the estimated matrix flow and preferential flow end member values for temperature, whereas one-sample Wilcoxon tests were used for a similar comparison of EC end members. Although all events showed that baseflow EC was significantly different from the EC in surface runoff, temperature end members were not significantly different in matrix flow and preferential flow for all events. Only events with distinct EC and temperature end members were used in further analysis. Durbin-Watson tests for autocorrelation of estimates of preferential flow were calculated for both methods using the lmtest package in R (Durbin & Watson, 1971;Zeileis & Hothorn, 2002).
Additionally, ion concentrations were compared between tile baseflow and surface event flow using mixedeffects models with repeated measures for events. Similar significance values were obtained using Mann-Whitney tests of event means. Most individual ions did not have significantly different end members; however, these were still compared against EC. In order to determine the ions most responsible for driving EC, individual mixed-effects linear models were used to determine the correlation between EC and tile ion concentration during events using the nlme package in R (Pinheiro, Bates, DebRoy, & Sarkar, 2018). In these models, EC was modeled against ion concentrations and site, with event as a random categorical variable. This allowed for within-event errors to be correlated, potentially accounting for uncontrolled differences between events due to event drivers (i.e., melt vs. storm), antecedent conditions, or management practices before the event (e.g., tilling or fertilizer application). Additionally, major ions, weighted by the magnitude, were summed into a single metric based on Equation D described by Visconti, dePaz, and Rubio (2010). Statistical analyses for this charge balance comparing end members and modeling against EC were conducted as described for individual ions.
Since known values of preferential flow could not be used for direct comparison, this study did not attempt to calibrate estimates of preferential flow or model it. Instead, our goal was to compare two simple options that could be applied widely across soil types, seasons, and event types. Estimates of preferential flow quantity at hourly and event scales were compared between the temperature and EC methods. An ANOVA was used to analyze the effects of soil (i.e., clay and loam), hydrograph position (i.e., rising and falling limbs), and season (i.e., growing season May-September and nongrowing season). A Breusch-Pagan test was also used to analyze heteroscedasticity of the estimates between comparisons. Estimates of preferential flow using both methods were also compared with TP and SRP loads at both hourly and flow-weighted event scales using ANOVA. A statistical threshold of α = .05 was used for all tests.

Event dataset
The observed events accounted for ∼58% of annual flow at the clay site, and 38% of annual flow at the loam site during the study period . The events covered a wide range of event types (e.g., snowmelt without rain, convection storms, rain on snow), event sizes, and seasons and were therefore an accurate and strong representation of climate and processes in the Great Lakes region. Figure 2 shows the rainfall, surface and tile flow, surface and tile temperature, tile EC, and estimated portion of preferential flow for six representative events. Supplemental Tables S1 and S2 provide more information for each event in the study for the clay and loam sites, respectively.

Temperature method
Sufficient temperature measurements were recorded in surface and tile flow for 38 of the events. However, when comparing tile temperatures during baseflow periods (i.e., 12 h before an event and 24 h after an event) with surface temperatures during the event, only 31 of the events showed significant differences between the two assumed end members (Supplemental Figure S1, Supplemental Tables S1 and S2). More of the events in the loam site had significantly different end members compared with the clay site (92 vs. 77%). During all events, surface flow ceased before tile flow returned to baseflow conditions, which necessitated the preferential flow end member for temperature to be estimated. Although the preferential flow end member estimate can be used to calibrate the preferential flow estimation method, there are likely many factors that can influence the end member value, making the end member estimate difficult and highly influential. Temperatures in event surface runoff ranged more widely than in baseflow or event flow in tiles. This suggests that estimating a constant end member would be F I G U R E 2 (a) Rainfall and (b) flow during six representative events at the two study sites separated by source (tile flow as solid lines and surface flow as dashed lines) and hydrograph position (baseflow in orange, rising limb in green, and falling limb in purple). Events were named by the site, then two-digit year, month, and day. Data for estimating preferential flow using (c) temperature and (d) electrical conductivity (EC) include tile (solid line) and surface (dashed line) measurements of temperature and EC, as well as the estimated values of matrix flow (dotted line). Estimated percentages of (e) preferential flow by both the temperature and EC methods are shown (blue and pink, respectively). Preferential estimates for temperature were noisier and values occasionally ranged outside of end members. However, the methods generally showed similar patterns within events a poor representative of preferential flow, and instead a dynamic end member value is necessary to accurately predict preferential flow. Similarly, temperatures ranged widely throughout the study from 0 to 23 • C in tile flow and from −1 to 21 • C in surface flow. During cold periods, tile flow was typically warmer than surface flow (e.g., Figure 2, Event 5), and the reverse was true during hot periods (e.g., Figure 2, Event 4). Most of the events without significant differences in temperature end members occurred during spring or mild winters when soil and air temperatures were more similar. This was similar to the findings of O' Driscoll and DeWalle (2006), where difference between end mem-bers, in this case air and stream water temperature, were greatest during winter events. During three events, some tile temperatures were outside of the temperature range between baseflow before and after the event (matrix flow end member) and surface flow during the event (preferential flow end member). However, all events had a least one time point when the event tile flow temperature was outside of the end member ranges for the given time point. This range was narrower than the range of baseflow temperatures because of the linear interpolation. With flow from events often lasting for several days, these out-ofrange values only accounted for 16% of temperature data. This likely coincided with periods later in events when the contribution of preferential flow is likely small (Vidon & Cuadro, 2010).
The dominant pattern in temperature within events was diel cycling of the surface temperature. This often resulted in surface temperatures oscillating about the tile temperatures (e.g., Figure 2, Event 4), even for events when the temperature end members were statistically distinct similar to oscillations in hyporheic temperature identified by Arrigoni et al. (2008). This caused estimates of preferential flow to oscillate as well, jumping to greater than one (i.e., when surface flow temperature dropped between tile flow and the estimate of matrix flow) and then to less than zero (i.e., when surface flow dropped below the estimate of matrix flow). An example of this is shown in Figure 2c, Event 4. However, the noisy and diel nature of temperature made it easier to visually distinguish the impact of surface flow on tile flow in the clay site. For events with generally high preferential flow, the tile temperatures mirrored the surface temperatures even after surface flow had stopped. This connection was less apparent in the loam soil where estimates of preferential flow were much lower. Overall, end member analysis using temperature estimates for matrix and preferential flow seemed to allow for a rough estimate of preferential flow quantity with a moderate number of exceptions due to indistinct end members and noise within events.
It is worth noting that temperature is not a conservative tracer, and convection of heat to groundwater and conduction to the soil likely further obscured the results (O'Driscoll & DeWalle, 2006). Indeed, thermal energy from surface runoff entering the soil could have moved directly into the tile water (advection) or been absorbed by the soil (convection). As the ratio of these two heat transfers is given by the Péclet number, estimating Péclet numbers for different sites and antecedent conditions may suggest when temperature delineation of preferential flow would be effective. Based on calculations at the field scale as described by Lyon and Troch (2007), Péclet numbers for the clay and loam sites were 0.42 and 22.2, respectively. However, this analysis considers heat transfer at the hillslope scale and not the soil profile that preferential flow encounters. Assuming macropores function as pipes, channeling thermal energy through the soil, we also estimated ranges of the Péclet number based on macropore velocities described by Gao et al. (2018) and macropore diameters observed by Grant, Macrae, and Ali (2019). The range in Péclet number for the clay site was 1.0 to 192, whereas the loam site ranged from 0.8 to 147. These wide ranges reflect the wide range in soil conditions during the study period (wet and dry, cold and hot) and make comparison between sites difficult.

Electrical conductivity method
Electrical conductivity measurements were recorded in surface and tile flow for 28 of the events. Observed EC of matrix flow during baseflow conditions (i.e., tile drainage) was significantly greater than the EC of the assumed preferential flow end member (i.e., surface runoff) for all events. However, since EC observed in baseflow was consistent and very different from EC observed in surface grab samples, we can be particularly confident in this method's distinct end members. The consistency even between sites suggests that this method could be easier to implement on a wide scale without much site-specific calibration being required. Tile EC responded quickly to event flow, in some events dropping before the hydrograph delineation method indicated that the events had begun (Figure 2, Events 2, 4, and 6). The responsiveness, especially in the clay soil, aligns with literature on rapid connectivity caused by macropores in dry soil (Graham & Lin, 2011). Conductivity even responded in the absence of surface runoff during snowmelt and behaved similarly to events that did have runoff (e.g., Figure 2, Event 5). The ability for this method to be applied to smaller events that do not generate surface runoff allows for a wide applicability of EC as a naturally occurring tracer for preferential flow separation. However, the slow return to pre-event EC, as seen by the steadily rising EC during baseflow after events (e.g., Figure 2, Events 1 and 3), suggests that surface flow is mixing with the matrix and lowering its EC. Cano-Paoli et al. (2019) found a similarly slow return to pre-event EC values in streams. On average, it took ∼3.5 d for the EC to return within a standard deviation of pre-event levels at the clay site and ∼0.5 d to recover at the loam site. Some events also showed estimated preferential flow at 100% persisting for hours at the hydrograph peak. Although this was likely event flow, it is unlikely that this represented the preferential flow we were trying to delineate. Indeed, in Pellerin et al. (2008), the highest percentage of new water (similar to preferential flow in our study) was 97%, with most events peaking around 60%. Higher predictions in our study may be attributed to the relatively high EC for the preferential flow end member (100 µS m −1 ) compared with values around 25 µS m −1 used by other studies (such as Pellerin et al., 2008). Although other studies use EC values from rainfall, we used EC values from surface runoff. Determining the appropriate EC value for preferential flow, as well as how much EC of matrix flow is affected by an event, likely requires more nuance and calibration than the simple linear interpolation used here.
Additionally, the smoothness of EC during events would suggest a consistent driver. However, when compared against individual anion and cation concentrations, EC showed generally weak linear correlations (Supplemental Figure S2). These correlations were much weaker than those observed by Chikhaoui et al. (2008), although their study did not contain as wide a distribution of event types and sizes. Ammonium and SO 4 2− had negative correlations with EC, and the correlation between K + and EC was near zero. These results are similar to the Michaud et al. (2018) study, which was conducted in similar climate and soils to this study. When considering the myriad sources and processing of individual ions, it is more understandable why correlations were weak (Laudon & Slaymaker, 1997). Ammonium, NO 3 − , K + , SO 4 2− , and Cl − were all applied at some point during the study at each site, though not always together. Additionally, NO 3 − , NH 4 + , and SO 4 2− availability change with biological activity (Russow, Stange, & Neue, 2009). Other ions, including CO 3 2− , Ca + , and Mg 2+ , are geogenic and are released by the calcareous soils in southern Ontario. Finally, anaerobic cycling affects the charges of some ions, such as Fe 2+ (not measured in this study) and CO 3 2− , which then affect EC differently. When event was included as a random effect in a mixed-effects linear model, correlations improved substantially, showing that individual ions do correlate with EC within events but not between events (Supplemental Figure S2).
The charge balance of major ions had a high correlation both without (r 2 = .64) and with (r 2 = .83) events included in the linear model (Figure 3, Supplemental Figure S2). This suggests that, although individual ions have poor relationships with EC, combining the information of multiple ions together does relate with EC and also explains the smoothness of EC both during and between events. Indeed, a linear model including all ions explained 83% of the variance in EC, but a Type II ANOVA only showed two statistically significant coefficients (SO 4 2− and K + ). End members were not significantly different for most ions or the charge balance, but this is likely due to high variability in sample concentrations due to processes and management. Given that the specific constituents driving EC are unknown, estimates of preferential flow using EC as a tracer should be treated with some caution. However, overall, EC seemed to provide a strong estimate of preferential flow that would be widely applicable.

Comparison of methods
Estimates of preferential flow generally agreed between the method using temperature and the method using EC (Figures 2e and 4). The plots of estimated preferential flow percentage during events were similar in shape for the majority of events, and the ranges also were typically F I G U R E 3 Electrical conductivity (EC) plotted against charge balance, a sum of major ions weighted by the magnitude of ion charge, measured in tile flow during eight events at the two study sites (clay in red and loam in blue). The correlation between EC and charge balance was r 2 = .64 based on a linear model and increased to r 2 = .91 when event was included as a random effect in a mixed-effects linear model similar. This suggests that the driver of the two methods is the same, likely preferential flow, and therefore both methods provide adequate approximations of actual preferential flow. However, when estimates based on EC and temperature during events were compared visually, the noisy nature of the method using temperature was more pronounced. Indeed, Durbin-Watson statistics for each method showed greater autocorrelation for the EC method (0.09) compared with the method using temperature (0.34), and both metrics were significant and suggest positive autocorrelation (i.e., preferential flow estimates at a timestep are correlated with the estimates at the previous timestep). Based on hourly data, there were very strong correlations in estimates of preferential flow rate using the two methods at the clay site, although the relationship was weaker at the loam site (Figure 4a). Correlation coefficients between hourly estimates of preferential flow using the two methods were .83 and .31 for the clay and loam sites, respectively. Similar results were found at the event scale (Figure 4b), where the loam site had weaker correlations between preferential flow estimates using temperature and EC when compared with the clay site (.89 and .98, respectively). As discussed above, EC was often more responsive early in events but took longer to return to pre-event levels, whereas tile event flow temperature oscillated more widely Vadose Zone Journal F I G U R E 4 Comparison of estimates of preferential flow based on temperature and electrical conductivity (EC) (a) at an hourly step (flow rate of preferential flow as q pref ) and (b) at the event scale (flow volume of preferential flow as Q pref ) at the clay and loam study sites (in red and blue, respectively). Samples within events are separated by the rising and falling limb of the hydrograph (in solid line with filled dot and dashed line with open dot, respectively). Events are grouped by season with the growing season (GS) shown in squares with tight dashes and the nongrowing season (NGS) shown in triangles with dotted lines. Flows are plotted on logarithmic axes, which emphasize lower flows and dampen higher flows, and the black diagonal line shows 1:1 estimates. Although there were generally strong correlations between estimates by each method, correlations at the clay site were stronger at both hourly and event scales. The rising limb and growing season at the loam site both had negative correlations between preferential flow estimates using the two different methods after surface flow ceased. Contrary to expectations based on these observations, the correlations between preferential flow estimates using the two methods were not markedly different on the rising and falling limbs of the event hydrographs at the clay site. However, there was a difference between the correlations on the rising and falling limbs at the loam site. The falling limb had a stronger correlation and the rising limb had a negative correlation between preferential flow estimates using the two different methods. This could be due to the wider range of estimates and larger sample size for the falling limb, which was always longer than the associated rising limb. However, the negative correlation in the rising limb suggests that one method may be erring and not just imprecise at estimating preferential flow on the rising limb. It is also possible that we did not have sufficient temporal resolution to accurately capture the rapid changes in preferential flow on the rising limb. By comparison, the falling limb would have had slower changes in preferential flow that were better captured by the 15-min data collection interval.
The clay site showed almost no impact of seasonality with correlation coefficients of .98 and .99 for preferential flow estimates using the two methods in the growing and nongrowing seasons, respectively. Although the data for the loam site may suggest some seasonality, there were not sufficient events within the growing season with both temperature and EC estimates of preferential flow to draw any definitive conclusions.
Both hourly and event-scale data showed weaker correlation at lower flows. Indeed, the residuals compared with a 1:1 line show significant heteroskedasticity, as confirmed by a Breusch-Pagan test. This suggests that aggregating across seasons or years will result in increasingly similar estimates of preferential flow due to the influence of larger events. Additional large events at the loam site may increase the weighting toward large events and strengthen the correlation between the two methods. Although a stronger correlation would result in more accurate total volume estimates of preferential flow, some important factors, such as soil wetting and hydrophobic interactions (Wang, Feyen, & Ritsema, 1998), may be drowned out during larger events and therefore are better understood during small events.

Preferential flow volumes
The clay site showed consistently higher percent preferential flow than the loam site ( Figure 5). However, within events, there were no significant differences between the rising limb and the falling limb at the clay site ( Figure 5a). Assuming preferential flow can be enhanced by desiccation cracks in clay soil, we would expect higher connectivity early in events that would decrease as the clay soil wetted (Grant, Macrae, & Ali, 2019;Hardie et al., 2012). It is possible that the temporal resolution of this study was not sufficient to accurately describe the rapid changes that may occur during the rising limb compared with the slower changes characteristic of the falling limb. In the loam soil where dry conditions are less associated with preferential flow, the difference between the rising and falling limbs was significant based on the temperature method, with the falling limb having more preferential flow (Figure 5b). This is similar to one of the sites described by Hardie et al. (2012). They posited this is due to restrictive deeper layers causing saturation at the soil surface, which causes connection of more macropores.  also found more preferential flow in loams under wetter conditions. Comparison between seasons showed no significant differences between the growing season (May-September) and the nongrowing season (October-April; Figures 5c and 5d). Plach et al. (2019) found that nongrowing season contributed significantly more to runoff than the growing season. Based on this, we would expect to see higher preferential flow during the non-growing season as well, but Vadose Zone Journal F I G U R E 5 Comparison of percentage preferential flow estimated using both temperature and electrical conductivity (EC) methods in the (a, c) clay and (b, d) loam sites. In Plots a and b, hourly data are split into rising limb (RL, purple) and falling limb (FL, green) of each event hydrograph. In Plots c and d, event data are split into the growing season (GS, May-September, blue) and the nongrowing season (NGS, October-April, orange). The yellow diamonds indicates mean of each subset. Only the comparison between RL and FL at the loam site was significant (indicated by *), with the falling limb of events showing a significantly greater percentage preferential flow than the paired rising limbs the results show no significant difference in percent preferential flow despite similar event sizes. In contrast, Grant, Macrae, Rezanezhad, and Lam (2019) found higher preferential flow in loam soils when frozen compared with unfrozen soil, although there were no differences in frozen and unfrozen clay soils. It is possible that milder winters in southern Ontario do not lead to frozen soils that facilitate preferential flow (Stadler, Stähli, Aeby, & Flühler, 2000), so that macropores function the same between the growing and nongrowing seasons. Taumer, Stoffregen, and Wessolek (2006) found lower effective cross section, a measure of area connecting flow through the soil, during summer and late fall, which may suggest that different processes associated with each season (i.e., drying and freezing) have similar effects on preferential flow.

Preferential flow and phosphorus load
Of the events analyzed for P loads, 11 events had sufficient TP and temperature data and 16 events had sufficient SRP F I G U R E 6 Comparison of estimated preferential flow and P concentrations during events at the clay and loam study sites (in red and blue, respectively). Phosphorus concentrations are separated into (a, b) soluble reactive P (SRP) and (c, d) total P (TP). Hourly P concentrations are compared against percentage preferential flow as estimated by (a, c) temperature and (b, d) preferential flow as estimated by electrical conductivity (EC). The data are subset by hydrograph position (rising limb in unfilled and falling limb in solid shapes) and by season (growing season [GS] in circles and nongrowing season [NGS] in triangles). The dashed blue line (b, d) highlights an observed threshold below which P loads at the loam site were often very low at ∼5% preferential flow and temperature for comparisons of temperature-based preferential flow predictions to P loads. Meanwhile, six and nine events had sufficient data for TP and SRP comparisons with EC estimates of preferential flow. Similarly to other studies comparing P and flow, both SRP and TP concentrations showed minimal correlation with preferential flow (Figure 6). Comparing loads and preferential flow volume resulted in stronger correlations, although these were still too weak to build predictive models (Table 1). There appeared to be a threshold at the loam site around 5% preferential flow, below which P concentrations were much lower (one to two orders of magnitude). This supports the hypothesis that P in tile drains is driven by connectivity with the surface. Additionally, stronger correlations of SRP compared with TP agrees with previous studies that show higher mobility of SRP through the soil (Jarvie et al., 2017;Michaud et al., 2018;Plach et al., 2019). This was only observed in the linear models and not the linear mixedeffects models (Table 1). TA B L E 1 Soluble reactive P (SRP) and total P (TP) loads compared with estimates of hourly preferential flow volume at two study sites. The R 2 values show the correlation between hourly P load and preferential flow volume based on a linear model (r 2 lm ) and a linear mixed-effects model with event as a random effect (r 2 lme ). Significance of the preferential flow terms of each model are shown as p lm and p lme for the linear model and the linear mixed-effects model, respectively Linear models between preferential flow and P loss that included event as a random effect resulted in correlations that were substantially higher, suggesting stronger relationships within events compared with between events (Table 1). Temperature-based estimates of preferential flow correlated more strongly with P than the EC method, which is likely due to the smaller number of events with sufficient data compared with the temperature method. Increased correlation when removing variation between events suggests that characteristics of the events themselves, or antecedent weather or management, impact P movement in addition to preferential flow or site (Vidon & Cuadra, 2010;Williams et al., 2016). Indeed, surface P concentrations ranged widely during the study, which likely obscured trends in comparing concentration and flow.
Also, as noted previously, the estimates of preferential flow in this study were not calibrated to measured preferential flow, which could explain this deviation. Underpredicting high P loads could be due to factors including recent fertilizer application or tillage (Williams et al., 2016) and leaching from frozen crops (Cober, Macrae, & Van Eerd, 2019), whereas overpredicting loads could be due to groundwater intrusion (Madison et al., 2014) or crops reducing soil disturbance (Ghestem et al., 2001). Erring more in larger events is likely to cause poor estimates of annual losses, so calibration efforts should emphasize these larger events. Due to the low number of events with sufficient data, we cannot draw strong conclusions at the event level but the generally positive relationship between preferential flow and P loss aligns with the hypothesis that preferential flow drives P loss. The weak correlation suggests that it is likely not the only driver. A possible driver in the calcareous soils at the study sites is supply limitation due to subsurface buffering by reactions with Fe 2+ (Plach et al., 2018), which was not measured in this study. This would likely not be a driver in soils saturated with P and also would be minimal due to flow in macropores interacting little with the subsoil. Further spatial and temporal event monitoring is necessary to average across all of this variability.

Extensions of this research
Based on the similarity between estimates of preferential flow using temperature and EC, further work is merited. Most field studies that measure tile flow also measure the temperature of the tile discharge; these data, coupled with measurements of surface water temperature or even air temperature, would allow for estimation of preferential flow across almost all tile-drained agricultural datasets. Although this study suggests these results are rougher than using EC, the volume of data would be sufficient for identifying major spatial and temporal trends in preferential flow across soil types. Coupling temperature-based estimates with EC estimates of preferential flow, which are less common, would allow for more fine calibration of estimates based on temperature. These estimates of preferential flow could be used to test or verify hypotheses posited by numerous studies about the relationships of preferential flow to soil structure , runoff volumes (Zhang et al., 2019), and P loads (Grant, Macrae, & Ali, 2019). This information could also be used to improve models of subsurface flow and P transport from soil columns (Urbina et al., 2019) to watershed and regional scales (Vogel, 2019).
As noted previously, this study did not directly measure preferential flow contribution to tile drainage because of the difficulties of applied tracers and large scales. The findings from this study confirm that temperature and EC both have potential for estimating preferential flow across soil types and seasons. Further work is necessary to calibrate these techniques from the current basic form to more accurately and precisely describe preferential flow. These should account for instances of poor estimation observed in this study as measured by disagreement between the methods, including diel temperature cycling, uncertainties about the ions driving EC, and nonlinear relationships. Other factors to consider include frozen ground (Stadler et al., 2000), saturation (Graham & Lin, 2011), no-till management (Williams et al., 2016), and plants (Ghestem et al., 2001). In estimating P transport due to preferential flow, reactions with and absorbance to soils should also be investigated.

CONCLUSIONS
This study compared estimates of preferential flow contribution with total tile drainage across several years in two soil types in southern Ontario. Both techniques were able to estimate preferential flow for most events, and the similarity in estimates suggests that these techniques could be effective for large-scale, widespread application in tile-drained agriculture. Estimates of preferential flow also correlated with P loads within events, which is further evidence that preferential flow contributes disproportionately more to P export from tile drained agricultural sites. Although further work is necessary to develop and calibrate these methods, this study shows that these rough techniques can be applied to data already collected to model and test posited explanations of observed phenomena in P, other nutrients, and water transport from tile-drained agricultural land.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.