Tile Drainage Increases Total Runoff and Phosphorus Export During Wet Years in the Western Lake Erie Basin

Artificial subsurface (tile) drainage is used in many agricultural areas where soils have naturally poor drainage to increase crop yield and field trafficability. Studies at the field scale indicate that tile drains disproportionately export large soluble reactive phosphorus (SRP) and nitrate loads to downstream waterbodies relative to other surface and subsurface runoff pathways, but knowledge gaps remain understanding the impact of tile drainage to nutrient export at watershed scales. The Western Lake Erie Basin is susceptible to summertime eutrophic conditions driven by non-point source nutrient pollution due to a shallow mean water depth and land use dominated by agriculture. The purpose of this study is to analyze the impact of tile drainage on downstream discharge, nutrient concentrations, and nutrient loads for 16 watersheds that drain to the Western Lake Erie Basin. Daily discharge and nutrient concentrations were summarized annually and during the main nutrient loading period (March–July) for 2 years representing normal nutrient loading period precipitation (2018) and above normal precipitation (2019). Results indicate positive correlations between watershed tile drainage percentage and runoff metrics during 2019, but no relationship during 2018. Additionally, SRP concentration and load were positively correlated to watershed tile drainage percentage in 2019, but not in 2018. Watershed tile drainage percentage was correlated with nitrate concentration and load for both years. The SRP concentration-discharge relationships suggested relatively weak, chemodynamic behavior, implying a slight enriching effect where SRP concentrations were greater at higher stream discharge conditions during both years. In contrast, nitrate concentration-discharge relationships suggested strong, enriching chemodynamic behavior during 2018, but chemostatic behavior during 2019. The difference in SRP and nitrate export patterns in the 2 years analyzed highlights the importance of implementing appropriate best management practices that target specific nutrients and treat primary delivery pathways to effectively improve downstream aquatic health conditions.

Artificial subsurface (tile) drainage is used in many agricultural areas where soils have naturally poor drainage to increase crop yield and field trafficability. Studies at the field scale indicate that tile drains disproportionately export large soluble reactive phosphorus (SRP) and nitrate loads to downstream waterbodies relative to other surface and subsurface runoff pathways, but knowledge gaps remain understanding the impact of tile drainage to nutrient export at watershed scales. The Western Lake Erie Basin is susceptible to summertime eutrophic conditions driven by non-point source nutrient pollution due to a shallow mean water depth and land use dominated by agriculture. The purpose of this study is to analyze the impact of tile drainage on downstream discharge, nutrient concentrations, and nutrient loads for 16 watersheds that drain to the Western Lake Erie Basin. Daily discharge and nutrient concentrations were summarized annually and during the main nutrient loading period (March-July) for 2 years representing normal nutrient loading period precipitation (2018) and above normal precipitation (2019). Results indicate positive correlations between watershed tile drainage percentage and runoff metrics during 2019, but no relationship during 2018. Additionally, SRP concentration and load were positively correlated to watershed tile drainage percentage in 2019, but not in 2018. Watershed tile drainage percentage was correlated with nitrate concentration and load for both years. The SRP concentration-discharge relationships suggested relatively weak, chemodynamic behavior, implying a slight enriching effect where SRP concentrations were greater at higher stream discharge conditions during both years. In contrast, nitrate concentration-discharge relationships suggested strong, enriching chemodynamic behavior during 2018, but chemostatic behavior during 2019. The difference in SRP and nitrate export patterns in the 2 years analyzed highlights the importance of implementing appropriate best management practices that target specific nutrients and treat primary delivery pathways to effectively improve downstream aquatic health conditions.

INTRODUCTION
Increased severity and frequency of Harmful Algal Blooms (HABs) caused by elevated toxinproducing phytoplankton and cyanobacterial biomass occurred in a majority of lakes around the globe from 1984 to 2012 (Ho et al., 2019). While HABs can be caused by various anthropogenic activities on the land, e.g., wastewater treatment plant discharges, septic system leaks, fertilizer runoff from residential lawns, fertilizer or manure runoff from farm fields and runoff or leaks from animal agriculture manure lagoons are the main HAB drivers in agricultural landscapes. Ultimately, severe HABs are costly for both human and ecological systems and can negatively impact drinking water sources, wildlife, tourism, recreation, property values, and commercial fishing. Annually, freshwater HABs are estimated to result in $4.6 billion of economic loss in the U.S. alone (Kudela et al., 2015). Direct HAB remediation is costly and typically involves either biological, physical, or chemical control measures; thus, mitigating nutrient losses at the source in agricultural fields could prove to be a cost-effective alternative to reduce HABs in receiving waters and lakes. It seems that a combination of both hydrologic and nutrient management practices are likely needed to reduce nutrient losses from agricultural fields and improve downstream aquatic health (Hanrahan et al., 2019).
Excess nutrients and HABs are not just problems for small inland lakes but can have significant impacts for large lake systems such as the Laurentian Great Lakes in the United States. Consider Lake Erie, the southern-most, shallowest, and most productive of the Laurentian Great Lakes, which contains <2% of the water held in the Great Lakes, but produces 50% of all fish harvested from the Great Lakes (Hushak et al., 1988). Summer-long HABs in Lake Erie can result in $5.6 million in lost fishing revenue, the world's largest walleye fishery (Wolf et al., 2017). The Western Lake Erie Basin (WLEB) supports about 80% of the sportfishing activity (Hushak et al., 1988) and is the most susceptible of the three Lake Erie basins to eutrophic conditions. The emergence of HABs in the WLEB can be partially attributed to the shallow depth (mean depth 7.3 m; ODNR, 2018) resulting in warm summer water temperatures, conducive to phytoplankton and cyanobacterial growth, and a drainage basin land use dominated by agriculture (71% of WLEB land area; USDA, 2005), which often transport excessive nutrient loads to the WLEB, particularly during the main nutrient-loading period occurring between March and July (e.g., Stumpf et al., 2012;Baker et al., 2019).
The most severe HABs occurring in the WLEB since 2002, indicated by the Lake Erie Harmful Algal Bloom severity index, have occurred in wetter than average years. Previous studies have identified very strong correlations between March-July stream discharge and March-July SRP and particulate phosphorus loads (R 2 = 0.83 and 0.90, respectively) from the major tributaries of the WLEB from 2002 to 2019 (Guo et al., 2020), implying phosphorus export is transport limited as indicated by chemostatic behavior with stream discharge (i.e., lower variability in concentration compared to stream discharge) (Williams et al., 2016a). From 1975 to 2017 annual precipitation in the Maumee River Basin, the largest tributary of the WLEB, increased by 102 mm, with increases in heavy (25.4-76.2 mm) and very heavy (>76.2 mm) precipitation events accounting for a majority of the increase, resulting in increased stream discharge and runoff ratio (defined as the ratio of area-normalized stream discharge to precipitation) in a majority of monitored streams (Williams and King, 2020). If this pattern of increased total rainfall and rainfall intensity continues, severe HABs will likely continue to be an economic and ecological problem for the region.
Recent studies have concluded that changes in SRP concentration vs. discharge (C-Q) relationships in WLEB tributaries, reflecting higher concentrations for a particular discharge compared to historic values, have been responsible for about two-thirds of the increase in SRP loads into the WLEB since the early 2000s, while the remaining one-third can be attributed to increased stream discharge (Jarvie et al., 2017;Choquette et al., 2019). Interestingly, changes in C-Q relationships for total phosphorus and three nitrogen species (total nitrogen, nitrate-nitrogen, and total Kjeldahl nitrogen) in WLEB tributaries were much less pronounced or showed decreases in C-Q relationships, resulting in little or no change in total load for these nutrient species (Choquette et al., 2019). Increases in SRP C-Q relationships have been hypothesized to be related to unintended consequences of long-term, large-scale changes in land management involving reduced tillage to decrease losses of particulate phosphorus, and increased artificial subsurface (tile) drainage in order to increase field trafficability and crop production while minimizing surface runoff (Jarvie et al., 2017). Collectively, these management practices have been suggested to have cumulative and converging impacts that have contributed to increased SRP loads to the WLEB.
Even today, tile drainage continues to be installed throughout the Midwest and the WLEB with many agricultural watersheds estimated to have over half their land under tile drainage. The large-scale prevalence of tile drainage has been shown to significantly alter the hydrology, and subsequent nutrient transport, at field and watershed scales. In Iowa, for example, tile drainage was shown to homogenize streamflow metrics across both high and low discharge conditions (Boland-Brien et al., 2014). The authors found that heavily tile-drained watersheds did not show the relationships between watershed size and streamflow metrics expected from natural systems (e.g., decrease in peak streamflow and increase in baseflow proportion with watershed drainage area). Other studies from Iowa indicated that tile drainage primarily affects the baseflow portion of the hydrograph (Schilling and Helmers, 2008a) and water exiting tile drains is primarily sourced from baseflow (Schilling and Jones, 2019). Contrasting results were recently reported in Ohio watersheds, whereby watersheds with extensive tile drainage had significantly lower streamflow recession constants and baseflow percentages, suggesting a flashier streamflow behavior in watersheds with high percentages of tile drainage (Miller and Lyon, 2021). The opposite pattern observed in Ohio watersheds compared to Iowa was hypothesized to be due to increased precipitation, decreased drainage capacity of soils derived from clay-and silt-sized glaciolacustrine deposits, and shallower water tables in Ohio compared to Iowa. Because of these differences in precipitation patterns and soil drainage properties, tile drains in Ohio tend to be installed at shallower depths and smaller lateral spacing (i.e., greater density) compared to installations in Iowa. In essence, tile drainage creates a more direct connection between precipitation, shallow subsurface water, and adjacent receiving streams in Ohio, compared to Iowa. What remains to be seen is the impact this more direct connection has on subsequent nutrient concentrations and loads to recipient freshwater systems -and the implications this could potentially have for management strategies to help alleviate HABs in the WLEB. As such, the purpose of this study is to evaluate the impact of tile drainage on watershed-scale hydrologic response and nutrient export behavior. Specifically, we determined the effect of tile drainage on (1) total runoff, (2) runoff ratio, (3) SRP and nitrate concentrations, (4) SRP and nitrate loads, and (5) C-Q dynamics annually and during the main nutrient loading period (March-July) for 16 watersheds located in the WLEB with varying proportions of tile drainage. For each watershed, tile drainage percentage was calculated from a recently developed 30m dataset (Valayamkunnath et al., 2020). Results were compared for two recent years, representing relatively normal (2018) and wet (2019) March-July conditions, which led to moderate, and severe algal bloom severity in the WLEB, respectively. Due to the strong surface connection that tile drainage creates (e.g., Smith et al., 2015;Jiang et al., 2021), we hypothesize that watersheds with the greatest percentages of tile drainage not only have greater SRP and nitrate concentrations and subsequent loads but also display more chemostatic C-Q behavior, especially during the main nutrient loading season (March-July).

Study Area
The WLEB is primarily contained within the Huron/Erie Lake Plain and Eastern Corn Belt Plains level III ecoregions (U.S. EPA, 2013). The Huron/Erie Lake Plain is a fertile, nearly flat plain with relic sand dunes, beach ridges, and end moraines. Historically, soil drainage was poor and elm-ash swamps, oak savannas, and beech forests were common. The 4,000 km 2 Great Black Swamp once existed in this ecoregion (Mitsch, 2017), but is now considered one of the most intensively artificially drained regions in the United States (Smith et al., 2008). The region now contains highly productive corn, soybean, and livestock farms due to the widespread installation of tile drainage. Bedrock topography is smooth in this ecoregion and is thought to have little impact on hydrology due to thick deposits of fine-grained clay-and siltsize glaciolacustrine sediments and glacial deposits of Wisconsin age (Ohio Division of Geological Survey, 2005; Macrae et al., 2021). The Eastern Corn Belt Plains are slightly hillier compared to the Huron/Erie Lake Plain and contain glacial deposits from Wisconsin-age end moraines. Soils in this ecoregion are loamier and slightly better drained compared to the Huron/Erie Lake Plain. Extensive corn, soybean, and livestock production replaced most of the original elm-ash swamps and beech forests.

Data
Daily stream discharge was downloaded from 16 United States Geological Survey (USGS) gages from 2018 to 2019 using the R package "dataRetrieval" (De Cicco and Hirsch, 2014) (Table 1; Figure 1). This time period was chosen to maximize the number of sites with available daily nutrient data in the WLEB and to align with the responses from the 2017 agricultural census data used to generate a tile drainage distribution map (Valayamkunnath et al., 2020). Daily streamflow was normalized by watershed area to calculate total daily runoff ("runoff "). Watershed characteristics and boundaries were obtained from the GAGES-II dataset (Falcone, 2011). Monthly precipitation data from 2018 to 2019 was aggregated to watershed boundaries to calculate monthly and annual precipitation (PRISM Climate Group, 2019).
These 2 years also had substantially different size of algal blooms develop in the WLEB as indicated by the Western Lake Erie Bloom Severity Index created by the U.S. National Oceanic and Atmospheric Administration's (NOAA) National Centers for Coastal Ocean Science (NCCOS). The severity index was created in 2002 and is based on the bloom's algae biomass over the peak 30-days. A goal of bloom severity ≤ 3 was established by NOAA and their research partners; however, the goal has only been met in 6 years since the index was created in 2002 and only once since 2008. For the years considered in this current analysis, the bloom severity index in 2018 was slightly higher than the goal (3.6), while the value in 2019 was much higher (7.5). An index above 5 indicates more severe blooms, and values > 7 are particularly severe, with extensive coverage affecting the lake. As such, this current effort targets a relatively low and high bloom severity year, respectively.
For each watershed, a 30-m resolution tile drainage map (AgTile-US) was used to calculate the percent of tile drainage occurrence (Valayamkunnath et al., 2020). This dataset used topographic slope, soil drainage, and county-level tile drainage census data to determine the most-likely tile-drained area in the contiguous United States with accuracy in the Midwest ranging from 82.7 to 93.6%. Cells with likely presence of tile drainage are Frontiers in Water | www.frontiersin.org indicated with a "1, " while no occurrence is indicated by a "0." The percentage of tile drainage was calculated by summing the total amount of tile-drained area by the total watershed area.
Daily nutrient concentrations (SRP and nitrate) were obtained from the National Center for Water Quality at Heidelberg University ("HU"; Tiffin, OH) and the USGS ( Table 1). Data from HU was downloaded from the Tributary Loading Website (https://ncwqr-data.org/HTLP/Portal), while the R package "dataRetrieval" was used to download daily USGS nutrient data (De Cicco and Hirsch, 2014). All nutrient sites considered are collocated at stream gages to allow for the calculation of nutrient loads. When multiple concentrations were reported for individual days, the mean, discharge-weighted, concentration was calculated.

Precipitation and Runoff Metrics
Monthly precipitation and daily runoff were summarized for each watershed annually, during the main nutrient loading period (Mar-Jul), and the months following the main nutrient loading period (Aug-Dec) for both 2018 and 2019. Runoff ratios were calculated by dividing runoff by precipitation, expressed as a percent. In addition, mean daily runoff for the highest 5% (Q05) and the lowest 5% (Q95) of annual observations was calculated for each watershed in 2018 and 2019. Annual and March-July runoff, and runoff ratios, were compared to watershed percentage tile drainage using Pearson's correlation. We used the Wilcoxon signed-rank non-parametric hypothesis test to evaluate differences in runoff, nutrient concentrations, and nutrient loads between 2018 and 2019, since the Shapiro-Wilk determined that the data were not normally distributed.

Nutrient Concentrations and Loads
Mean discharge-weighted SRP and nitrate concentrations were calculated annually, during the main nutrient loading period (Mar-Jul), and the months following the main nutrient loading period (Aug-Dec) for each year. Daily nutrient loads were calculated by multiplying stream discharge by the nutrient concentration and summarized annually, during the main nutrient loading period (Mar-Jul), and the months following the main nutrient loading period (Aug-Dec).

Concentration-Discharge Characterization
The transfer of nutrients and water through a watershed is inherently related. Nutrient concentrations at watershed outlets reflect differences in application rates, flow paths, residence times and can be considered a "fingerprint" of watershed transport, reaction, and mixing processes (Knapp et al., 2020). Calculating the variation in nutrient concentration as function of stream discharge therefore provides insight into how watersheds store and release water and nutrients. We quantified SRP and nitrate C-Q relationships two separate ways: (1) by fitting power-law relationships between nutrient concentrations and discharge (e.g., Godsey et al., 2009;Musolff et al., 2015;Knapp et al., 2020) and (2) by calculating the ratio of coefficient of variation (CV) for nutrient concentrations and discharge to assess the relative variability of each (Thompson et al., 2011;Gorski and Zimmer, 2021). Nutrient concentration-discharge (C-Q) behavior was estimated by fitting the following power-law relationship: which is identical to the linear relationship in doublelogarithmic space: log 10 (C) = log 10 (a) + b * log 10 (Q) where log 10 (a) and b are the intercept and slope of the C-Q relationship, respectively (Godsey et al., 2009;Knapp et al., 2020). Chemostatic C-Q behavior is defined as b ≈ 0 and implies little to no variation of C with changes in Q, while a value of 1 or −1 suggests nutrient concentrations vary proportionately positively or negatively to discharge, respectively. Nutrients with b = 1 display higher concentrations at high discharge conditions (i.e., mobilization), while nutrients with b = −1 indicate lower concentrations at high discharge conditions (i.e., dilution). However, b ≈ 0 does not indicate that the variability in concentration is small and the goodness of fit of the power-law C-Q relationship is uninformative as b approaches zero (Thompson et al., 2011). The slope of the power-law C-Q relationship (b) does not reveal whether changes in concentration, discharge, or a combination of both were responsible for the observed relationship. For these reasons, an additional, non-parametric measure of chemostatic behavior was calculated defined by the ratio of the coefficients of variation (CV) between concentration and discharge: where µ and σ represent the mean and standard deviation, respectively. Chemostatic C-Q behavior is observed when CV C is much smaller than CV Q resulting in CV C /CV Q << 1. Both C-Q metrics (b and CV C /CV Q ) are consistent with each other, but CV C /CV Q is more general and does not imply nutrient concentrations are truly invariant, or uncorrelated with changes in discharge (Thompson et al., 2011).

Watershed Characteristics
Mean watershed area (Figure 1) was 3,550 km 2 and ranged from 11 to 14,248 km 2 ( Table 1). The average percent of watershed area under tile drainage was 45.3% and ranged from 30.2 to 62.6%. Agriculture, which includes pasture and cultivated crops, was the dominant land use for all watersheds and ranged from 67 to 86% (mean = 78%).

Precipitation and Runoff
Annual and March-July precipitation and runoff ( Table 2) were significantly greater in 2019 compared to 2018. Both years had significantly greater annual precipitation and runoff compared to 1991-2019 averages; however, March-July precipitation and runoff were similar in 2018 compared to 1991-2019 averages, while March-July precipitation and runoff were significantly   Table 2).
Tile drainage was significantly positively correlated to annual runoff (Pearson's r = 0.52), March-July runoff (Pearson's r = 0.70), annual runoff ratio (Pearson's r = 0.53), and March-July runoff ratio (Pearson's r = 0.60) in 2019 but was not related to any runoff metrics in 2018 (Figure 3). Neither August-December runoff or runoff ratio was found to be correlated with tile drainage.

Nutrient Concentrations and Loads
SRP and nitrate concentrations and loads are summarized for the 16 study watersheds annually, during the main nutrient loading period (Mar-Jul), and the months following the main nutrient loading period (Aug-Dec) for 2018 and 2019 in Table 3. Mean annual, March-July, and August-December concentrations of SRP and nitrate were similar between the 2 years. Mean annual and March-July nitrate loads and annual SRP loads were similar between the 2 years; however, there was a significantly greater March-July SRP load in 2019 compared to 2018. The SRP and nitrate loads from August-December were significantly greater in 2018 compared to 2019.

Concentration-Discharge Characterization
The slopes of the SRP concentration-discharge relationship were similar during the 2 years when comparing annual and March-July data (Table 4). However, the slopes of the nitrate C-Q relationship were significantly greater in 2018 compared to 2019.
Annual and March-July SRP coefficient of variation ratios were similar for SRP. The annual coefficient of variation ratios for nitrate were similar for the 2 years, but the March-July coefficient of variation ratio was significantly greater in 2018 compared to 2019.
Watershed tile drainage (%) was not related to the slope of the March-July SRP or nitrate C-Q relationship in 2018 or 2019 ( Figure 5). However, the SRP coefficient of variation ratios were significantly negatively correlated to tile drainage using annual (Pearson's r = −0.65) and March-July (Pearson's r = −0.72) data in 2019, but not in 2018. The nitrate coefficients of variation ratios were not related to tile drainage in either year using both annual and March-July data.

Connections Between Tile Drainage, Stream Discharge, and Nutrient Transport
In 2019, which had above average annual precipitation and March-July precipitation, watershed tile drainage percentage was significantly correlated with annual and March-July runoff and annual and March-July runoff ratios (Figure 3). These results are supported by recent research from Miller and Lyon (2021) who observed a flashier streamflow response in watersheds with high percentages of tile drainage in Ohio compared to watersheds with little to no tile drainage. The lack of relationship between tile drainage and annual or March-July runoff and runoff ratios during 2018 suggests that the presence of tile drainage amplifies the flashiness during wetter than average years in this region, while the effect is less noticeable in years with normal (or low) amounts of precipitation. Given the trend in more frequent large rain events occurring in the WLEB (Williams and King, 2020), we suspect that tile drainage will likely contribute to flashier downstream behavior in the future. Management efforts that reduce the flashiness of rivers following storms events such as drainage water management, drainage water recycling, vegetated buffers, or wetland restoration are likely to reduce the amount of nutrient-rich water leaving farm fields and improve downstream aquatic health.
There was a noticeable difference between March-July SRP and nitrate concentrations and loads between for the 2 years analyzed. The wetter than average conditions during March-July of 2019 resulted in significant correlations between tile drainage and SRP mean concentrations and loads (Figures 4A,B). However, during 2018 there were no relationships between tile drainage and mean SRP concentrations and loads (Figures 4A,B). This suggests that tile drainage exacerbates the transport of SRP in particularly wet conditions. In fact, at the field-scale, numerous studies have documented how tile drainage is connected to the surface through macropores and other preferential flow paths (Smith et al., 2015;Williams et al., 2016a;Macrae et al., 2019). Storm events have been found to accelerate the transport of particulate and dissolved nutrient species in tile drainage (Jiang et al., 2021). Significantly greater precipitation and runoff occurred during the months following the main nutrient loading period (Aug-Dec) in 2018 compared to 2019. Tile drainage was found to be significantly correlated with SRP and nitrate concentrations and loads during this period in 2018, but not in 2019. This supports evidence that watersheds with large percentages tile drainage accelerate the transport of nutrients during wet conditions, but do not amplify nutrient export during drier conditions. Repetitive applications of fertilizer or manure to soils has led to an accumulation of phosphorus at the soil surface and pronounced degrees of vertical soil phosphorus stratification whereby soil phosphorus levels decrease greatly with depth (Robbins and Voss, 1991;Sharpley, 2003). Vertical soil phosphorus stratification is generally more noticeable in no-till or conservation-till soils compared to conventional tilled soils (Vu et al., 2009;Sharpley et al., 2012). Initial recommended fertilizer levels were established based on a long history of conventional tillage, whereby phosphorus would be expected to be incorporated into the plow-layer (0-20 cm). However, currently, a majority of fertilizer is not incorporated into the plow layer, but broadcast onto the surface, which has likely contributed to increases in SRP into the WLEB in recent decades (Smith et al., 2015. This underscores the importance for measuring soil test phosphorus concentrations to identify fields that are at risk for greater phosphorus loss (Duncan et al., 2017). During wetter than average years, shallow groundwater tables are elevated relative to dryer years and near-surface macropores are activated more often resulting in more potential phosphorus loss from heavily stratified soils (Williams et al., 2016b).
In both 2018 and 2019 tile drainage was significantly correlated with March-July nitrate concentrations and loads, despite significantly greater precipitation occurring in 2019 (Figures 4C,D). The strong correlation between nitrate concentrations and loads with tile drainage suggests that watersheds with greater percentages of tile drainage export larger nitrate loads regardless of precipitation conditions. Nitrate concentrations are generally thought to be greater in groundwater compared to surface runoff. Schilling and Helmers (2008b) found that 70% of nitrate in tile-drained fields was transported via diffuse flow, while suspended sediment and phosphorus transport were more associated with faster flow regimes. This helps explain the similar nitrate concentrations recorded in 2018 and 2019, despite more rain occurring in 2019.
Tile drainage is designed to reduce surface runoff and erosion. While particulate phosphorus and erosion have likely reduced  partially due to widespread installation of tile drainage, the significant increase in dissolved phosphorus in the last 20 years has been proposed as an unintended consequence of an increase in tile drainage installations and reduced tillage designed to improve soil stability (Jarvie et al., 2017). In fact, research suggests that tilling soil may be beneficial in breaking up macropores that are more likely to develop under reduced or no-till conditions (Williams et al., 2016a).  We explored the relationship between various watershed characteristics, runoff, and nutrient concentrations to determine how variation in watershed properties affected nutrient transport. Interestingly, agricultural land use (percentage of watershed area) was not found to be correlated with annual or March-July SRP concentrations for either year, but significantly positively correlated with nitrate both annually and during March-July.
Widespread accumulation of nitrogen in the root zone of agricultural soils from decades of fertilizer application is likely a major reason why tile drainage and agricultural land were highly correlated with nitrate in both 2018 and 2019 (Van Meter et al., 2016). Whereas, the transport of SRP is confined to shallower soil layers that are activated only in wetter years.

Variability of Nutrient Concentration-Discharge Dynamics With Tile Drainage Intensity
Annual and March-July SRP displayed relatively weak, enriching chemodynamic behavior as indicated by low positive SRP C-Q slopes in both 2018 and 2019 (b = 0.26-0.34; Table 4; Gorski and Zimmer, 2021). The lack of difference in SRP C-Q slopes between the 2 years, both annually and during March-July, suggests that transport of SRP is more strongly influenced by precipitation and stream discharge than C-Q dynamics. In fact, mean March-July SRP loads were significantly different between the 2 years, but not mean March-July concentrations (Table 3).
Interestingly, annual and March-July nitrate C-Q dynamics displayed different behavior in 2018 and 2019 (Table 4), despite similar mean nitrate concentrations and loads between the 2 years ( Table 3). Annual and March-July nitrate C-Q slopes had stronger, enriching chemodynamic behavior in 2018 (b = 0.72). However, the annual nitrate C-Q slope in 2019 showed a moderate, enriching chemodynamic behavior (b = 0.24), while March-July nitrate C-Q slopes indicated chemostatic behavior (b = 0.01).
Chemostatic nitrate C-Q behavior has previously been demonstrated in systems dominated by agriculture (e.g., Basu et al., 2010;Thompson et al., 2011), due to consistent nitrate availability from contributing areas that have uniform nitrate concentrations. However, as the contributing area shrinks following precipitation events, near-stream areas with greater denitrification rates and lower nitrate concentrations contribute relatively more to stream discharge compared to upland areas with greater nitrate concentrations (Marinos et al., 2020), resulting in contributing areas with more heterogeneous nitrate concentrations and a more chemodynamic nitrate C-Q behavior. For examples Marinos et al. (2020) and Gorski and Zimmer (2021) found chemodynamic nitrate C-Q behavior to be dominant at low flows, while chemostatic nitrate C-Q behavior was more common during high flow events. These observations are consistent with our results showing more chemodynamic nitrate C-Q behavior in 2018, compared to the chemostatic response observed in 2019.
All annual and March-July coefficient of variation ratios for SRP and nitrate displayed moderate, enriching chemodynamic behavior as indicated by values > 0.3 (Thompson et al., 2011). The mean annual calculated coefficient of variation ratio for SRP from the 16 watersheds was very similar to the results obtained by Williams et al. (2016a) from two watersheds in the WLEB (0.52 and 0.54). The coefficient of variation ratios for SRP and nitrate aligned with most of the SRP and nitrate C-Q slopes, except for March-July nitrate C-Q behavior in 2019, which indicated chemostatic behavior according to the C-Q slope. This lack of agreement between the two metrics suggests that there are other reasons for concentration variation beyond variability in stream discharge. The coefficient of variation ratio does not explain the variance in concentration, but avoids assuming a particular underlying concentration-discharge relationship (Thompson et al., 2011).
The general lack of relationship between watershed tile drainage (%) and nutrient C-Q parameters eludes to other management factors effecting solute concentration dynamics. Marinos et al. (2020) found a significant correlation between tile drainage and nitrate C-Q behavior in the Upper Mississippi Basin, where watersheds with high percentages of tile drainage spent a greater amount of time displaying chemostatic behavior. However, the distribution of tile drainage in the watersheds analyzed in their study was much greater than those compared in this study and included many sites with little to no tile drainage. The watersheds compared in this study had at least 30% of their land area underlaid by tile drainage and thus the true effect of tile drainage on nutrient C-Q dynamics is perhaps not attainable with the current selection of sites without a comparison to watersheds with little to no tile drainage. Unfortunately, most of this region has been extensively tiled (Figure 1) to increase agricultural yields and thus unimpacted watersheds with water quality and discharge records do not exist.

Potential Limitations and Concluding Remarks
Clearly, data availability may impact the interpretation of our results since only 2 years were compared. However, including more years of stream data would have resulted in many fewer watersheds to compare. In addition, the tile drainage dataset used to create watershed-scale estimates of tile drainage coverage is based on data collected in 2017, thus including older stream data in our analysis may not be appropriate. Further, watersheds with lower percentages of tile drainage were only monitored more recently meaning these sites would be excluded if the analysis were extended further back in time. Regardless of the limited amount of data available, the 2 years chosen for analyses offer stark differences in nutrient loading period precipitation and subsequent WLEB algal bloom severity. From this perspective, these 2 years offer contrasting conditions and insights when we, for example, consider catchment management practices that target nutrient load reductions.
Specifically, there was significantly greater precipitation in the WLEB during 2019 for the primary nutrient loading period compared to 2018, which experienced a normal level of precipitation during this time. During this wetter year we identified significant positive correlations between watershed tile drainage density, runoff and SRP concentrations that were not found in the drier year. From a phosphorous management perspective, this difference in response suggests that tile drainage may amplify nutrient export during years with wetter than average precipitation. When considering how we manage for nitrogen at the watershed-scale, tile drainage density was significantly correlated with nitrate in both years. Such constant importance of drainage density implies precipitation amounts may be less important for nitrate transport during this time. Further, the concentration-discharge metrics were, in general, not correlated with tile drainage density and suggested moderate, enriching, chemodynamic behavior. Nitrate concentration-discharge metrics appeared to be more affected by the different precipitation amounts in 2018 and 2019, whereas SRP concentration-discharge metrics were similar in 2018 and 2019. Collectively, these relationships would point to nitrogen being more impacted by deeper legacy sources compared to SRP being related or confined to more shallow soil layer sources. Understanding these dynamic relationships between hydrologic response and nutrient concentrations in streams forms a cornerstone to adequate design and management capable of limiting loading across a range of climatic forcing. This latter aspect is increasingly important as weather extremes and climatic variability increase globally.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
SM designed the study, performed the analyses, and prepared the first draft of the manuscript. SL secured funding for the project, advised the interpretation of the data, and helped draft and edit the manuscript. Both authors contributed to the article and approved the submitted version.