Tile drainage as a driver of streamflow flashiness in agricultural areas of the Midwest, USA

The Midwest of the USA is a highly productive agricultural region, in part due to the installation of perforated subsurface pipes, known as tile drains that remove excess water from wet soils. Tile drains rapidly move water to nearby streams, influencing how quickly streamflow rises and falls (i.e., streamflow “flashiness”). Currently, there are no comprehensive studies that compare the extent to which tile drainage influences flashiness across large and diverse agricultural regions. We address this knowledge gap by examining growing‐season (April–October) flashiness using the Richards‐Baker Index (RBI) in 139 watersheds located throughout the Midwest. Using a spatial tile‐drainage dataset, watersheds were split into low, medium, and high tile‐drainage classes. We found no significant differences between the flashiness of these three classes using a one‐way Kruskal–Wallis test. When watersheds were separated into infiltration groups to help control for different soil types, the high tile‐drainage class RBI was significantly higher than the low tile‐drainage class RBI in the high infiltration group. To further understand the causes of flashiness, additional environmental variables and their relationship to flashiness were examined using multivariate regression. In the low infiltration group, tile drainage significantly reduced flashiness, with watershed area and average depth to water table being the largest influences on flashiness. Tile drainage produced a larger reduction in flashiness in the high infiltration watersheds, with the largest influences being percent clay in the watershed and watershed area. These results indicate that the influence of tile drainage on flashiness emerges only after other watershed variables are accounted for. Given that tile drainage may increase in the future as precipitation patterns and extremes change, flashiness will likely continue to be modified. These results lead to an improved understanding of flood‐generating and nutrient transport mechanisms that are relevant to stakeholders across a wide range of sectors.


| INTRODUCTION
As one of the most productive agricultural regions in the world, the Midwest of the USA produced over 63 billion dollars worth of crops in 2018 (USDA National Agricultural Statistics Services, 2017) in areas converted from the natural wetland, prairie, and forest ecosystems to agricultural lands (Radeloff et al., 2005;Slater & Villarini, 2017).Historically, agriculture in many parts of the Midwest was limited by high water tables.The installation of agricultural tile drainage removed excess water in these areas, allowing for cultivation to occur (Blann et al., 2009;Spaling & Smit, 1995).Tile drainage is an agricultural practice that places perforated pipes under the surface of a field at depths of 0.6-1.2m to lower the water table and quickly remove excess soil water from the field, resulting in drier soil that promotes larger rooting zones, earlier planting dates, and higher yields (Blann et al., 2009;Schilling & Helmers, 2008;Spaling & Smit, 1995).The installation of tile drains transformed the way that watersheds respond to precipitation inputs (Raymond et al., 2008;Slater & Villarini, 2017) and, due to their ability to rapidly transport nutrients, lowered the water quality in many bodies of water.The impacts extend to areas as large as Lake Erie and the Gulf of Mexico, resulting in ongoing eutrophication and hypoxia (Blann et al., 2009;Goolsby et al., 2001;McIsaac & Hu, 2004;Randall & Mulla, 2001).
Tile drains lead to altered streamflow regimes due to subsurface flow that causes water to move from soil to stream faster than in unmodified conditions (Blann et al., 2009;Rodvang & Simpkins, 2001;Schilling et al., 2015).Annual baseflow, which is the portion of streamflow that comes from groundwater and sustains river discharge during dry periods, has been found to increase due to tile drainage in Iowa (Boland-Brien et al., 2014;Schilling, 2005;Schilling & Helmers, 2008;Schilling & Libra, 2003), while tile drainage decreased baseflow in Ohio (Miller & Lyon, 2021).Additionally, tile drainage has been found to contribute 30%-61% of annual streamflow in small headwaters in Ohio (King et al., 2014), while in Iowa tile drainage was found to contribute 15%-43% of annual streamflow (Arenas Amado et al., 2017).
In the Minnesota River Basin, Foufoula-Georgiou et al. (2015) found that tile drainage has increased streamflow at the intermediate percentiles, suggesting tile drainage may affect typical (normal) streamflow more than streamflow extremes.
As a result, changes in flashiness can occur due to alterations in precipitation or other watershed characteristics, especially at the land surface, that modify the partitioning of precipitation into surface runoff or infiltration, both of which control how quickly water reaches a stream channel (Gannon et al., 2022;Liu et al., 2006;Zhang & Schilling, 2006).
While the influence of tile drainage on streamflow flashiness has been studied at the watershed or state level, how tile drainage influences flashiness at the regional spatial scale has not been examined.
We found only two studies that examine the influence of tile drainage on flashiness in multiple watersheds.The first (Boland-Brien et al., 2014) uncovered that watersheds (n = 24) in Iowa without tile drainage were flashier than those with tile drainage, possibly due to topography.Tile drains are usually installed in flatter, poorly drained areas, so watersheds with less tile drainage tend to have larger slopes that can drain water more rapidly, creating flashier conditions.Additionally, they found that tile-drained watersheds with high infiltration soils exhibited increased baseflow and reduced flashiness (Boland-Brien et al., 2014).The second study, performed in Ohio, found that sites (n = 59) with more tile drainage tended to be flashier than sites with less tile drainage (Miller & Lyon, 2021).These conflicting results indicate that tile drainage may not always lead to streamflow flashiness, with changes perhaps related to local conditions such as topography and soil type (Miller & Lyon, 2021).Additionally, Gannon et al. (2022) found that flashiness was influenced by subsurface characteristics across much of the Midwest, but this study only posited that the effect of subsurface characteristics on flashiness was due to tile drainage because of its prevalence in the region.
Iowa and Ohio generally represent the east-west limits of the highly tile drained regions of the Midwest; however, there is little information about tile drainage and its influence on flashiness between these two states.This work fills a knowledge gap in understanding the relationship between flashiness and tile drainage in Midwestern watersheds.Our primary research question is, after controlling for other watershed characteristics, how does tile drainage influence flashiness at the scale of the Midwest?While high flow events are essential for aquatic and riparian ecosystems (Baker et al., 2004;Poff et al., 1997), further understanding flashiness and its drivers is necessary.Streams that are flashier can experience increased flooding (Wiskow & van der Ploeg, 2003) and a decline in water quality, both of which are harmful to the environment and economy (Miller & Lyon, 2021;Rabalais et al., 2002;Rabotyagov et al., 2014).

| Site selection
Watershed boundaries, stream-gaging stations, and daily streamflow data from 2010 to 2019 were retrieved from the United States Geological Survey (USGS) National Water Information System (NWIS; USGS, 2016) for the Midwest (defined in this study as Minnesota, Iowa, Wisconsin, Illinois, Indiana, Michigan, and Ohio; Figure 1).This time period was chosen to align with that of the tile drainage dataset (additional details in Section 2.2).The streamflow data used for the analysis were limited to the warm season (April-October) to reduce the influence of rain-on-snow or frozen ground precipitation events, which can lead to large runoff and streamflow pulses (Lapenta et al., 1995) that could result in streamflow flashiness not related to tile drains.Watershed characteristics such as LULC data, the number of dams and their distance from the gages, soil type, slope, and climate characteristics were extracted from GAGES-II (Falcone, 2011) and the 2016 National Land Cover Dataset (NLCD; Dewitz, 2019) for each watershed.NLCD data was used rather than land cover data from GAGES-II because it was more recently updated and matched the time period of the tile-drainage dataset (see Section 2.2).Following Miller and Lyon (2021), only watersheds that were >25% agricultural and <25% urban were included (Dewitz, 2019).Additionally, USGS gages had to have fewer than six major dams (dams ≥15 meters in height) in the watershed, have a watershed area of less than 2000 km 2 , and be >8 km downstream from the nearest major dam (Falcone, 2011), and have <10% missing data (Miller & Lyon, 2021).These criteria reduce the possibility of other factors, such as dams or urban areas, from obfuscating the influence of tile drainage on flashiness in the study watersheds.Based on these criteria, 139 USGS gaging sites were selected across the Midwest (Figure 1; Table S1).

| Flashiness metric
The (Matlab) Toolbox for Streamflow Signatures in Hydrology (TOSSH; Gnann et al., 2021) was used to calculate the Richards-Baker Index (RBI; Equation 1).RBI was developed as a flashiness metric and has been applied in the Midwest (Baker et al., 2004) and elsewhere in the United States (Dow, 2007;Gannon et al., 2022;Jarvie et al., 2017;Saxe et al., 2018): In Equation 1, q i is streamflow on day i, q iÀ1 is streamflow on the previous day, and n is the total number of days.Given that RBI is the sum of the absolute differences of daily streamflow on two successive days divided by the sum of the total streamflow for April to October (Baker et al., 2004), it is a dimensionless (contrast) metric that is standardized by the flow magnitude over the sample period.Flashier streams, therefore, have a higher RBI due to larger daily differences than in streams with more stable flow.

| Watershed characteristics
Spatial tile-drainage data for this study was gathered from AgTile-US, a 30 m resolution raster dataset that indicates the presence or absence of tile drains in each given cell across the continental United States (Valayamkunnath et al., 2020; Figure 1).This dataset was created using a combination of United States Department of Agriculture (USDA) Tile Census data, National Land Cover Database (NLCD) cropland, slope, and Soil Survey Geographic Database (SSURGO) soil drainage categories.Across the Midwest, it was found to have an accuracy range of 83%-94% (Valayamkunnath et al., 2020).The number of tile-drained cells within each watershed were summed and used to calculate the percent area of tile drainage (Miller & Lyon, 2021).Following this, the watersheds were separated into three classes based on the percent area of tile drainage: low (5-15% drainage; n = 22), medium (>15%-40% drainage; n = 51), and high (>40% drainage; n = 66).Binning by tile drainage rather than working with continuous percent area is useful for broad categorical comparisons and also helps to control for the uncertainty and inaccuracies (e.g., tile drain density, spatially limited validation with satellite data) in the AgTile dataset.However, tile data was treated as continuous in the multivariate regression (discussed below) in order to quantify the influence that it has on flashiness.
Additional watershed characteristics were retrieved from the GAGES-II dataset.We first extracted the percent of the watershed that had soils in hydrologic groups of A, B, C, or D from GAGES-II, which classify soils based on infiltration potential (Falcone, 2011).
Watersheds with a majority of type A or B soils were separated into the "high" infiltration group (n = 68), while watersheds that had a majority of C or D soils were separated into the "low" infiltration group (n = 71).Watersheds with soil types in these groups have differing responses to precipitation, which, if analysed together, can dampen the flashiness signal (Boland-Brien et al., 2014;Miller & Lyon, 2021).Other GAGES-II data that may also influence flashiness included watershed area (km 2 ), mean watershed slope (%), % of the watershed with clay soils present, and depth to water table (seasonally high; m).Gridded (4-km resolution) growing season ETₒ (mm), represented as a short reference crop calculated with the ASCE Penman-Monteith Method (Walter et al., 2005), and precipitation (mm) data were retrieved from GridMET (Abatzoglou, 2013) for 2010-2019.The climate data were summed into a seasonal total (April-October) from the daily data, which were then averaged at the watershed level to gain a seasonal average of precipitation and ETₒ for each watershed.Alongside tile drainage, several of these data (area, slope, clay percentage, depth to water table, and climate characteristics) have been found to influence flashiness (Boland-Brien et al., 2014;Miller & Lyon, 2021).

| Statistical analyses
The effect of tile drainage on streamflow on the three (low, medium, and high) tile-drainage classes was examined using a Kruskal-Wallis test.Type-I error of α = 0.05 was used to evaluate statistical significance for this and all following statistical analyses.Spearman's correlation coefficient (r s ) was used to test for correlations between tile drainage and flashiness at sites in Ohio and Iowa, while a two-tailed Kolmogorov-Smirnov (K-S) test was used to compare the distributions of RBI for the three tile-drainage classes with and without dams.
A multivariate linear regression was used to further understand the connections between climate, watershed characteristics, and RBI.
Multivariate regressions with variables such as ETₒ, precipitation and land use regressed on streamflow have been widely used to model complex hydrologic processes due to the various factors that influence streamflow (Ficklin et al., 2018;Jiang et al., 2011;Thomas, 2000;Zhang et al., 2016).In this study, RBI is the response variable while multiple watershed and climate variables and tile drainage are the predictor variables.Predictor variables were z-score standardized to control for differences in their magnitudes (Sirdas & Sen, 2003).Standardized regression coefficients also allow us to compare the relative effect size of different predictors.The accuracy of the model was evaluated with the coefficient of determination (R 2 ), mean absolute error (MAE), and the distribution of residuals.MAE was used as a metric of average error due to it being a more interpretable measure of mean error (Willmott & Matsuura, 2005) while providing an error metric in the units of the response variable.near or outside of glacial boundaries (Naylor et al., 2016).While most of the flashiest watersheds were in the medium and high tile-drainage classes, there was no spatial pattern of these classes having consistently higher RBI (Figure 2), although many of the higher values are present in areas that had more tile drainage in the eastern portion of the study area (see Figure 1).Many of the medium and high tiledrainage watersheds, especially in the northern and western portions of the Midwest (Figure 2), exhibited low RBIs that were similar to those in the low tile-drainage classes.Additionally, most of the low drainage sites had RBI that were similar to those of nearby medium and high drainage watersheds, suggesting that tile drainage is not a primary control on flashiness at this scale.
There was little difference in the RBI distributions (Figure 3), medians, and means of the low (n = 22), medium (n = 51), and high (n = 66) tile-drainage classes, respectively.The tile-drainage classes had similar RBI distributions and, using the Kruskal-Wallis test, were found to have no significant differences in their distributions (p = 0.36), suggesting no clear connection between tile drainage and flashiness at these sites.

| RBI distributions by soil infiltration group
In order to better determine how tile drainage affects flashiness, the sites were divided into low and high infiltration soil groups (Figure 4a).
Using the Kruskal-Wallis test, we found no significant difference (p = 0.18) between the RBI distributions of the three tile-drainage classes in the low infiltration drainage group.However, the high tile- drainage class exhibited a median value that was higher than that of the medium and low tile-drainage classes, indicating that some influence of tile drainage may have been present.
In contrast to the watersheds with low infiltration soils, the sites with a majority of high infiltration soils were less flashy (Figure 4b),   1;

| Multivariate regression
Figure S2), with a visual inspection of the residuals indicating that they were approximately normally distributed for both models.
In contrast to the findings from the tile-drainage class analysis, the multivariate regression found that tile drainage had a negative influence on flashiness in both infiltration groups when accounting for other variables (Table 1).The magnitude of the effect of tile drainage on RBI in the high infiltration soils was approximately twice the strength of its influence in low infiltration soils.In the low infiltration soils, depth to water table and watershed area were the two strongest influencing factors, while precipitation was the third strongest factor.
For high infiltration watersheds, the largest regression coefficient was the percentage of clay in the watersheds.The second largest influence on RBI was tile drainage, which, as with the low infiltration sites, had a negative influence on RBI (i.e., tile drainage reduces flashiness).
Alongside clay and tile drainage, watershed area, slope, and precipitation were all influencing factors on tile drainage at these sites.
Between the two groups, precipitation and watershed area had the strongest influence on the low infiltration group due to the runoffgenerating capabilities of those soils.In the high infiltration soils, the influence of clay was about three times greater than that in the low infiltration soils, while tile drainage and slope both had influences that were far larger than the low infiltration sites.

| Tile-drainage classes
No coherent spatial pattern of increasing flashiness based on tile drainage class was found, indicating that local factors might be exerting more influence on flashiness across the Midwest.Some of the areas with the highest RBI, despite having low tile drainage, were in the southern portions of the Midwest where there is a convergence between a lack of glaciation (Naylor et al., 2016) and low infiltration soils.These unglaciated regions tend to have shallower soils and steeper slopes in many areas, which are major factors in causing faster runoff responses and flashier streams (Baker et al., 2004).Agriculture within the unglaciated areas is concentrated in river valleys, which may allow for flashier behaviour due to the proximity of tile drainage to the stream channel, creating shorter travel times and a less attenuated precipitation response.In comparison to Miller and Lyon (2021), who included areas outside of glacial boundaries in their tile-drainage analyses, our study contains a smaller proportion of sites in these areas.Moreover, the multivariate regression models account for some of the factors (such as slope and soil type) that may be different outside of the glacial boundary and thus influence flashiness (Morton et al., 2015).
A large reason for a lack of difference between the tile-drainage classes was due to watersheds with majority low or high infiltration soils being analysed together.Hydrologic soil type can cause tile drained watersheds to respond to precipitation differently due to infiltration potential (Miller & Lyon, 2021), which had the effect of muting the influence of tile drainage on flashiness.When the watersheds were split into low and high infiltration groups (Figure 4), the results of the low infiltration group found no significant increases in flashiness between the medium and high tile classes, although the median value did increase, which aligned with findings from Miller and Lyon (2021) for Ohio.However, the high infiltration group diverged from the findings in Boland-Brien et al. ( 2014), who found that flashiness T A B L E 1 Results of the multivariate linear regression, including the number of sites, regression performance metrics, and (standardized) coefficients for the variables.decreased with increasing tile drainage, as flashiness significantly increased from the low to high tile-drainage class.This is likely due to our use of sites with >5% tile drainage, as sites with <5% tile drainage have been found to have higher slopes and, therefore, an influence on flashiness (Boland-Brien et al., 2014).Removal of sites with <5% tile drainage likely reduced the flashiness of the low tile-drainage class, which therefore allows flashiness to increase with tile drainage as some of the influence of slope is removed.Additionally, many sites in the high infiltration group had higher percentages of low infiltration soils than the Iowa sites from the Boland-Brien et al. ( 2014) study, which may have contributed to tile drainage increasing flashiness.
Another factor that may have influenced the results in these analyses is area.As the area of a watershed increases, it generally becomes less flashy.The low tile-drainage classes had the smallest median area for all sites (Figure 3) and high infiltration sites (Table S2; Figure 4b), which allows for some bias to be present when comparing the tiledrainage classes.

| Influences on flashiness
Major dams are an additional factor present in some watersheds, although our analysis limited them in both their number and distance from the gaging station.Dams are known to alter the hydrograph by reducing peak flows and increasing low flows, which can reduce flashiness (Baker et al., 2004;Graf, 2006).Although our original methods controlled for the number and distance of major dams from the gages, we performed an additional analysis on sites without major dams (n = 116) to test for an influence due to damming.The results with major dams removed from the low (n = 18), medium (n = 40), and high (n = 58) tile-drainage classes showed small differences in the flashiness distributions, with no significant differences found ( p = 0.54; p = 0.42; p = 0.82) when compared to the original tiledrainage classes with the two-tailed K-S test.The differences between the without-dams median compared to the median from the original results were 0.108, 0.014, and 0.049 for the low, medium, and high tile-drainage classes, respectively.This, combined with the lack of difference between the distribution, indicated that there was little effect on the median RBI of these watersheds due to major dams.
The true amount of tile drainage in each watershed, compared to the modelled amount from AgTile (Valayamkunnath et al., 2020), may also have influenced the results.While the AgTile dataset was ground-truthed with 16 000 points at sites across the Midwest (Valayamkunnath et al., 2020), these sites were scattered and may not be representative of all the tile drained areas in this region.The AgTile dataset also does not estimate the density of tile drains within each 30-m grid cell.Tile drainage systems with more closely spaced drains can quickly transport more water from the soil to nearby streams or rivers (Kennedy et al., 2012;Schilling et al., 2015;Strock et al., 2010), so flashiness may have varied for watersheds with different tiledrainage density.Sites with shallower tile drains could also increase flashiness by more directly connecting precipitation inputs with streamflow (Miller & Lyon, 2021), while sites with deeper installation depths would have slower response times to precipitation inputs and intercept more groundwater (Schilling et al., 2015), therefore reducing flashiness.While we are unable to confirm this from the AgTile dataset, the low infiltration sites in this study likely had tile drains installed closer to the surface due to shallower seasonally high water table depth (median: 0.7 m) when compared with the water table depth of high infiltration sites (median: 1.1 m; Falcone, 2011).Moreover, AgTile-US gives no indication of surface intakes or areas where controlled drainage is being implemented due to the difficulty in assessing where these practices are being used.Both practices could confound the results of this study by either reducing tile flow due to controlled drainage (Sunohara et al., 2015) or increasing tile flow due to surface inlets (Schilling & Helmers, 2008).

| Comparisons with previous work
The impact of tile drainage on streamflow flashiness has not been analysed for most watersheds within the Midwest, which limits comparisons across much of the region.Thus, we compared our flashiness results with those of the statewide Iowa (Boland-Brien et al., 2014) and Ohio (Miller & Lyon, 2021) studies, which represent the spatial ends of the highly tiled portion of the Midwest.The sites within these respective states were extracted and their RBI values were correlated with tile drainage.We found that RBI had a weak negative correlation to tile drainage in Iowa (r s = À0.13;p = 0.54; n = 26) and a moderate positive correlation to tile drainage in Ohio (r s = 0.53; p = 0.01, n = 24).In comparison, Boland-Brien et al. (2014) posited sites in Iowa to have decreased flashiness with increased tile drainage, as sites with tile drainage had lower mean flashiness.Our results in Iowa find a weak and insignificant negative correlation between tile drainage and flashiness, although the results cannot be directly compared because no correlation was calculated in the Iowa (Boland-Brien et al., 2014) study. In Ohio, Miller andLyon (2021) found that tile drainage was negatively correlated to the flashiness metric (r = À0.57;p < 0.05), which indicated that flashiness increased with tile drainage (although opposite sign from our correlation due to their use of a different flashiness metric: TQ mean ; Konrad & Booth, 2002).
Our correlation-based findings in the Ohio watersheds closely match those found in Miller and Lyon (2021).

| Connections between RBI and watershed variables
From the multivariate regression (Table 1), we found that tile drainage had the second-strongest relationship with RBI in high infiltration sites.In comparison, it had the fourth strongest relationship with RBI in low infiltration soils.Generally, the multivariate regression found that tile drainage has a negative influence on RBI in the high infiltration soil group, which aligns with previous studies of this soil type (Boland-Brien et al., 2014).The multivariate regression also found that tile drainage had a negative effect on flashiness in the low infiltration group, which contradicts previous work (Miller & Lyon, 2021) and our bivariate correlations in Ohio.This is likely due to the multivariate regression controlling for other variables that can influence flashiness (such as clay % and depth to water table) when assessing the role of tile drainage on RBI, which a bivariate correlation is unable to do.Moreover, the reduction in flashiness is likely a product of reduced soil moisture due to tile drainage (Lam et al., 2016), as well as an increase in baseflow in areas that have high infiltration soils (Boland-Brien et al., 2014).The Ohio sites in this study had a watershed average of 10.5% high infiltration soils, while all low infiltration sites have an average area of 15.7% high infiltration soils.This means that, especially near the boundary of where the high and low infiltration groups are, there is a possibility that mixed responses from the soil types are present, which may cause the signal of tile drainage to be altered.
Watershed area was included as a predictor variable in the multivariate regression due to the wide variety (7-1900 km 2 ) of watershed areas included in this study.As watershed area increases, flashiness tends to decrease due to larger watersheds having more baseflow, more varied characteristics, and variable inputs from tributaries (Baker et al., 2004;Gannon et al., 2022;Singh, 1997).Despite the influence that area has on flashiness, we did not limit the minimum watershed size because half of the sites that were less than 100 km 2 had flashiness values similar to larger watersheds.Additionally, the multivariate regression controlled for area when assessing the influence of other variables, such as tile drainage, on RBI.Watershed area had the second-strongest negative relationship with RBI in the low and fourth strongest in the high infiltration groups.Watershed area in the low infiltration group likely had a stronger negative influence on RBI due to the overall flashier nature of those watersheds, which would allow for a larger change in flashiness as watershed area increased and more varied inputs (streams and groundwater) attenuated flow.The sites had similar median areas of 527 km 2 for low infiltration watersheds, compared to a median value of 521 km 2 for the high infiltration watersheds, indicating that there was no large difference between their sizes.A characteristic that can be related to watershed area is slope, which was also included in the regression models.Smaller watersheds may also have larger slopes, although the multivariate regression found that slope had a negative influence on flashiness in both infiltration groups.In previous work, sites with little tile drainage have been found to have larger slopes and an increased influence on flashiness (Boland-Brien et al., 2014).Therefore, we removed sites that had less than 5% tile drainage to help control for slope.Slope may have had a negative influence on RBI since removing the sites with the least tile drainage did not remove all the sites with high slopes, and many sites with smaller slopes had higher flashiness due to other watershed characteristics.
The multivariate regression found the strongest positive relationship of all high infiltration variables was between clay percentage and RBI, while it had a small influence on RBI in low infiltration watersheds which was likely due to these watersheds having slightly higher percentages of clay present.The low infiltration sites had less variability in clay (%) than the high infiltration soils.The standard deviations and median clay percentages between tile classes were similar in the low infiltration soils, while the high infiltration sites had more variation in their clay percentages (Table S2).This may have allowed clay percentage to exert more influence on RBI than in the low infiltration sites, which had less variability.Moreover, the response of RBI likely varied between infiltration groups because the high infiltration sites had a higher baseline infiltration rate than sites in the low infiltration group (Dripps & Bradbury, 2010).This means that as clay increases in the high infiltration watersheds, there can be a larger infiltration reduction because of the larger difference between the infiltration rates of the high infiltration soils and clay.In areas that contain high infiltration soils, increased clay has been found to increase runoff due to the reduced ability of soil to infiltrate water (Anderson et al., 2009).
In contrast to clay, the strongest predictor of RBI in the low infiltration soils was depth to water table.While the high infiltration soils had deeper water tables (Table S2; Falcone, 2011), the low infiltration sites had a wider range of water table depths.However, these depths were on average shallower than the high infiltration soils, which likely increased the influence of depth to water table on RBI in the low infiltration soils.Sites with shallower water tables may not be able to infiltrate as much precipitation, increasing runoff and flashiness when compared to sites with deeper soils and water tables that can infiltrate more precipitation.Tile drainage is oftentimes installed in areas with shallow water tables (Blann et al., 2009), which creates a possibility that some effect from tile drainage may be connected to water table depth.This, combined with the low infiltration soils having shallower soils, creates the possibility that some effect from tile drainage may be hidden by the depth to water table.

| Connections between RBI and climate variables
Precipitation had the third strongest positive influence of the predictor variables in the low infiltration soils, while it had the fourth strongest influence on RBI in the high infiltration soils.Part of the reason for weaker relationships is due to the inherent variability of precipitation (and, perhaps, the use of seasonal total precipitation).Across the Midwest, the northern portion (Minnesota, Wisconsin) receives less precipitation during the growing season than areas that are more southern and eastern (Indiana, Ohio; Morton et al., 2015).This gradient was present in our study, with watersheds in southern Indiana and Illinois receiving over 800 mm of precipitation (April-October mean of the 10 years in our study), while watersheds in Minnesota, Wisconsin, and Michigan received mean seasonal total precipitation in the range of 570-700 mm over the course of this study.Combined with the variability in tile drainage, this precipitation gradient likely contributes to precipitation not having a stronger relationship with flashiness.Additionally, the flashiest sites in this study tend to be in the southern and eastern areas, which is likely due, at least partially, to the higher precipitation in those areas.The variability of the magnitude, intensity, and location of precipitation are all lost, especially in larger watersheds (Singh, 1997), since precipitation is a basin-wide and multi-year seasonal average in this study.Moreover, the amount of runoff generated by precipitation can also be influenced by antecedent soil moisture, which can increase the variability of an area's runoff response (Tarasova et al., 2018;Williams et al., 2023).The multivariate regression found that the strength of the relationship between ETₒ and RBI was weak and insignificant.While precipitation drives the rising limb of the hydrograph, ETₒ can be a controlling factor of the recession limb (Karlsen et al., 2019;Tashie et al., 2019).However, given that we used a seasonal average of ETₒ, much of this relationship is likely lost due to ETₒ remaining relatively constant despite the rising and falling of the hydrograph, whereas any change in precipitation will directly influence the rising limb of the hydrograph.
Despite its lack of influence on RBI (Table 1), ETₒ was not removed from the model because it did not cause the model to perform more poorly.Moreover, keeping ET o in the model allows for a more comprehensive understanding of factors that do and do not affect flashiness.

| Future directions
While we attempted to account for the main drivers of flashiness, other factors such as sinuosity, channel morphology, the density of tributaries, and riparian buffers could all influence flashiness in watersheds that otherwise have similar characteristics.For future studies, more inter-watershed variability could be controlled for with a geographically weighted regression (Brunsdon et al., 1996), which could account for the inherent spatial variability of climate, soil type, and other watershed characteristics that this study did not take into account.Moreover, further studying flashiness across wetter or drier seasons or years, rather than annual or growing-season averages, may help to better understand how antecedent soil moisture can influence runoff generation at a regional scale in a tile-drained landscape.Additionally, using data that is of higher temporal resolution may help to further understand the influence that tile drainage plays on streamflow flashiness and its change over time.

| CONCLUSIONS
Using 139 watersheds, we analysed the effect of tile drainage and other watershed and climate characteristics on flashiness in the Midwest.Tile-drainage data was gathered from the AgTile-US dataset and used to separate watersheds into low, medium, and high tile-drainage classes.Using a Kruskal-Wallis test on the three drainage classes, we found no influence of tile drainage present.When accounting for soil type, we found that there was a significant difference between the low and high tile-drainage classes in the high infiltration group, while there were no significant differences in the low infiltration group.
Given that many watershed and climate characteristics can influence flashiness, we developed a multivariate regression for the two infiltration groups.We found that, when controlling for other watershed variables, tile drainage reduced flashiness in both the high and low infiltration groups.This is likely due to tile drainage reducting soil moisture when the influence of other variables that affect flashiness is held constant.Other factors that are present in areas with tile drainage (such as % clay or shallower water tables) also influence flashiness, but their impacts are difficult to disentangle from tile drainage with univariate or bivariate analyses such as graphical summaries or correlation.However, uncertainties in the true amount of tile drainage likely reduced the accuracy of our Kruskal-Wallis test and regression results.While tile drainage reduced flashiness, our results indicate that other variables, such as % clay, depth to water table, slope, watershed area, and precipitation can act to mask the influence of tile drainage on flashiness.While many of these variables will stay constant in the future (i.e., % clay, slope, area), changes to precipitation patterns and extremes may drive an increase in installation (tiled area or tiled density) in order to transport larger quantities of water from agricultural areas.Such alterations to tile drainage would likely modify the relationship between flashiness and tile drainage further and contribute to water quality issues since flashiness and high flows are drivers of nutrient transport.Without further management efforts, such as controlled drainage, enlarged riparian buffers, and agricultural best management practices, precipitation changes and possible changes to tile drainage area and/or density could continue to exacerbate water quality issues across the Midwest.

1
United States Geological Survey (USGS) streamflow gages included in this study.As an example, USGS 04195500 (top right panel) is shown with the tile drainage dataset in the background.The tile-drainage class is portrayed using three symbols with circle for low drainage, square for medium, and triangle for high.Black symbols indicate watersheds with majority high infiltration soils, while red indicates watersheds with majority low infiltration soils.The 30-m resolution AgTile-U.S. dataset is shown in the background, with dark grey indicating areas that are estimated to have tile drainage.
| RBI distributions by tile drainage class RBI varied spatially across the 139 watersheds in the Midwest, with values ranging from 0.04 (low flashiness) to 1.08 (high flashiness) and a median value of 0.36 (Figure 2).RBI throughout the Midwest decreased with increasing latitude, where the higher latitude watersheds located in Minnesota, Wisconsin, northern Iowa, and Michigan exhibited the majority of the smaller (0-0.3)RBI values.The watersheds that exhibited the highest (0.6-1.08) RBI were located in Illinois, Indiana, and Ohio, where many of the highest values for low and medium classes were located in southern Illinois and Indiana in areas

F
I G U R E 2 Richards-Baker Index (RBI) for the study watersheds.Lighter symbols characterize watersheds with more stable hydrographs, while dark blue symbols portray watersheds with flashier hydrographs.The 30-m resolution AgTile-U.S. dataset is shown in the background, with dark grey indicating areas that are estimated to have tile drainage.
with tile-drainage class mean and median values well below those of the low infiltration tile drainage groups.For these high infiltration sites, the Kruskal-Wallis test found a significant difference (p = 0.038) between the RBI distributions of the low and high tiledrainage classes.With this group, the mean and median values generally increased as tile drainage increased, indicating that tile drainage may have had some effect on RBI, although the effect was not large enough for a statistically significant difference between the medium tile class and the other two classes.However, given that many watershed characteristics influence RBI, it is likely that there were other factors influencing flashiness besides tile drainage in both the low and high infiltration groups.

A
multivariate regression model was developed to better quantify the strength of the relationship between each watershed characteristic and flashiness while controlling for the other watershed characteristics.This model included RBI as the response variable, with six watershed characteristics (watershed area, slope, clay %, precipitation, ETₒ, and depth to water table) and tile drainage percentage as continuous predictor variables.The models resulted in high coefficients of determination in both the low (R 2 = 0.58; MAE = 0.105; p < 0.05) and high (R 2 = 0.65; MAE = 0.069; p < 0.05) infiltration groups (Table Boxplots comparing the Richards-Baker Index (RBI) distributions and means (red dots) of the low (n = 22), medium (n = 51), and high (n = 66) tile drainage.Grey dots are scaled to show watershed areas.The middle bar of the boxplot represents the median RBI, while the top and bottom of the box represent the 75th and 25th percentiles respectively.The end of the whiskers are the maximum and minimum values that are not considered outliers.F I G U R E 4 Boxplots of the tile drainage classes for low and high infiltration soils.The grey dots are scaled to show watershed areas, while the mean value for each class is shown by a red dot.The low infiltration tile drainage classes have 9, 31, and 31 sites, respectively.The high infiltration tile drainage classes have 13, 20, and 35 sites, respectively.The middle bar of the boxplot represents the median RBI, while the top and bottom of the box represent the 75th and 25th percentiles respectively.The ends of the whiskers are the maximum and minimum values that are not considered outliers.
values for the coefficients are given in parentheses.Statistically significant coefficients ( p < 0.05) are bolded.