The Role of Hydraulic Connectivity and Management on Soil Aggregate Size and Stability in the Clear Creek Watershed , Iowa

The role of tillage practices on soil aggregate properties has been mainly addressed at the pedon scale (i.e., soilscape scale) by treating landscape elements as disconnected. However, there is observed heterogeneity in aggregate properties along flowpaths, suggesting that landscape scale hydraulic processes are also important. This study examines this supposition using field, laboratory and modeling analysis to assess aggregate size and stability along flowpaths under different management conditions: (1) tillage-induced abrasion effects on aggregate size were evaluated with the dry mean weight diameter (DMWD); (2) raindrop impact effects were evaluated with small macroaggregate stability (SMAGGSTAB) using rainfall simulators; and (3) these aggregate proxies were studied in the context of connectivity through the excess bed shear stress (δ), quantified using a physically-based landscape model. DMWD and SMAGGSTAB decreased along the flowpaths for all managements, and a negative correspondence between the proxies and δ was observed. δ captured roughness effects on connectivity along the flowpaths: highest connectivity was noted for parallel-ridge-till flowpaths, where δ ranged from 0–8.2 Pa, and lowest connectivity for contour-ridge-till flowpaths, where δ ranged from 0–1.1 Pa. High tillage intensity likely led to an increase in aggregate susceptibility to hydraulic forcing, reflected in the higher gradients of aggregate size and stability trendlines with respect to δ. Finally, a linear relationship between DMWD and SMAGGSTAB was established.


Introduction
Agricultural management practices associated with tillage have been shown to impact soil aggregate size characteristics and stability through the mechanized weakening and breaking of the aggregate structure, i.e., organized groupings of different size soil grains and organic material with defined void spacing [1][2][3][4].The repeated disturbance to the soil increases the recovery time needed for the cohesion between soil particles and organic material to strengthen and build aggregate structure [5,6].During a season, the tillage disturbances can shift the fraction of the aggregate size distribution towards the lower size range of 0.25-2.0mm, hereafter referred to as the SMall macroAGGregate (SMAGG) [7] based on the hierarchical classification of Oades and Waters [5].Thus, differences in SMAGG size fractions have been found to be associated with differences in land use and management [8][9][10].
Tillage also reduces residue cover, leaving the soil surface exposed and vulnerable to rainfall and runoff processes [11,12].During a rainfall event, kinetic energy is transferred from falling raindrops into the soil aggregate structure upon impact [13][14][15].When the applied raindrop kinetic energy exceeds the internal molecular bond strength of the soil aggregate, the aggregate structure fails [16].Additionally, rainwater penetrates the void spaces of the aggregate structure, and, as the excess internal pore pressure builds up, it causes the aggregate to disaggregate [17].The disaggregated soil fragments can become more easily entrained by fluid forces and transported downslope [18,19].Entrained soil aggregates within overland flow are also further subject to hydraulic forces and breakdown during transport [20,21].
An extensive literature survey reveals that the following key approaches have been used to quantify the role of management on soil aggregate size characteristics and stability: (1) the use of dry sieving [22][23][24] in order to make inferences regarding the degree of abrasion to the soil triggered by tillage and, thus, its quality [25,26], and (2) the use of rainfall simulation or wet sieving tests, in order to account for the rate of wetting on the slaking and breakdown of soil aggregates [10,15,16,[27][28][29][30][31].Rainfall simulation tests are preferred against wet sieving tests, when trying to mimic breakdown under applied raindrop kinetic energy [6,17].Most of the aforementioned studies have focused on large macroaggregates or clod-sized fractions (>2 mm), not on the SMAGG fraction, even though the spatial redistribution of SMAGG fraction has been found to be most sensitive to management practices, such as tillage intensity, compared to larger size fractions [8][9][10]32].Therefore, there is a degree of uncertainty in their findings and their interpretation.
In addition, most of these aforementioned studies are limited to aggregate characterization at the soilscape scale, defined as a cluster of polypedons with homogenous composition of soil [19,33,34], by treating landscape elements as disconnected, i.e., independent from each other.Extension of their findings to actual landscapes, which consist of a series of interconnected soilscapes, is therefore challenging.Connectivity is influenced by management.The heterogeneity in landscape properties and fluxes controls the dynamic patterns of aggregate mobilization, transport, breakdown, and deposition by flow within the landscape [35][36][37][38][39][40][41][42][43].There are studies that have examined aggregate properties at different hillslope positions and found significant spatial heterogeneity in aggregate stability amongst positions-e.g., [44,45].However, these studies have not examined the effects of landscape hydraulic connectivity on aggregate stability along the same flow pathway with homogeneous management and soil texture, and, thus, the role of hydraulic connectivity on the variability in aggregate stability across a landscape is unknown.
The breakdown of aggregates in overland flow (in a laboratory flume) under different roughness conditions was studied by Wang et al. [21], who found that hydraulic properties controlling bed shear action, such as flow depth and friction factor, play a significant role in aggregate breakdown with distance.They also found that the degree to which the hydraulic properties affect breakdown is dependent on the aggregate size fraction.Their findings, along with the findings from the other studies mentioned above, suggest that hydraulic connectivity along with management practice dictate the size distribution and stability of aggregates found on the landscape.However, no study has systematically interlinked soil aggregate properties (size and stability) to flowpaths and hydraulic connectivity under different management practices.We posit that the degree of spatial heterogeneity in both aggregate size distribution and stability along a flow pathway is affected by both hydraulic connectivity and the type and intensity of tillage practices [11,12,[46][47][48][49] (tillage intensity herein refers to the number of tillage occurrences, and the type and severity of the disturbance i.e., tillage depth and percentage of disturbed area [50]-more frequent and severe disturbances to the soil are considered to be of high intensity, whereas less frequent and minimal disturbance is considered to be of low intensity).It is assumed therefore that the spatial heterogeneity of the aggregate proxies should reflect the spatial heterogeneity of hydraulic connectivity and management across dominant flowpaths of the landscape [45] (Figure 1).Addressing the above supposition requires a method for examining the role of hydraulic connectivity on aggregate size distribution and stability under different management practices.While rainfall simulators can be effectively used to capture the effects of raindrop-induced breakdown on SMAGG properties, they do not consider the role of hydraulic bed shear action along a flowpath and, thus, cannot be used to examine the role of hydraulic connectivity under different management practices [16].Physically-based landscape models [49,51,52] can be used to close this research gap by serving as tools to link the action of bed shear with SMAGG properties along a flowpath, provided that they can capture the feedback between flow and roughness (managementcontrolled) effectively.The enhanced WEPP model by Papanicolaou et al. [49] satisfies this condition and, thus, can be used to improve our understanding of SMAGG properties along a flow pathway under different management practices.
The goal of this study is to examine the size distribution and stability of SMAGGs at different hillslope locations and under different management practices, within the context of hydraulic connectivity.The study is performed at the Clear Creek watershed, a site of the Intensely Managed Landscapes Critical Zone Observatory in southeastern Iowa in the U.S. Midwest, where management practices related to tillage and surface roughness modification and their resultant effects on runoff allow us to study the roles of tillage and hydraulic connectivity on aggregate properties.The nature of the proposed study is unique, as it considers a holistic assessment of aggregate properties using field, laboratory and computational analysis: (i) we evaluate the role of tillage-induced abrasion on aggregate size distribution through measurement of the dry mean weight diameter (DMWD); (ii) we account for the raindrop impact kinetic energy through measurement of SMAGG stability (SMAGGSTAB) using rainfall simulation tests; (iii) finally, we consider the two aggregate proxies in the context of flowpath hydraulic connectivity through the bed shear stress, which is computed using the Papanicolaou et al. [49] modeling framework that encompasses landscape connectivity and roughness features (i.e., random roughness, oriented roughness, and row orientation) resulting from three different tillage regimes.Addressing the above supposition requires a method for examining the role of hydraulic connectivity on aggregate size distribution and stability under different management practices.While rainfall simulators can be effectively used to capture the effects of raindrop-induced breakdown on SMAGG properties, they do not consider the role of hydraulic bed shear action along a flowpath and, thus, cannot be used to examine the role of hydraulic connectivity under different management practices [16].Physically-based landscape models [49,51,52] can be used to close this research gap by serving as tools to link the action of bed shear with SMAGG properties along a flowpath, provided that they can capture the feedback between flow and roughness (management-controlled) effectively.The enhanced WEPP model by Papanicolaou et al. [49] satisfies this condition and, thus, can be used to improve our understanding of SMAGG properties along a flow pathway under different management practices.
The goal of this study is to examine the size distribution and stability of SMAGGs at different hillslope locations and under different management practices, within the context of hydraulic connectivity.The study is performed at the Clear Creek watershed, a site of the Intensely Managed Landscapes Critical Zone Observatory in southeastern Iowa in the U.S. Midwest, where management practices related to tillage and surface roughness modification and their resultant effects on runoff allow us to study the roles of tillage and hydraulic connectivity on aggregate properties.The nature of the proposed study is unique, as it considers a holistic assessment of aggregate properties using field, laboratory and computational analysis: (i) we evaluate the role of tillage-induced abrasion on aggregate size distribution through measurement of the dry mean weight diameter (DMWD); (ii) we account for the raindrop impact kinetic energy through measurement of SMAGG stability (SMAGG STAB ) using rainfall simulation tests; (iii) finally, we consider the two aggregate proxies in the context of flowpath hydraulic connectivity through the bed shear stress, which is computed using the Papanicolaou et al. [49] modeling framework that encompasses landscape connectivity and roughness features (i.e., random roughness, oriented roughness, and row orientation) resulting from three different tillage regimes.

Site Description
This study was conducted in southeastern Iowa within agricultural fields of the Clear Creek Watershed (CCW), Iowa (41.74 • N, −91.94 • W), which is part of the U.S. National Science Foundation Intensely Managed Landscapes Critical Zone Observatory [53] (Figure 2).In CCW, the mean annual precipitation is 889 mm/year and temperature is 9 • C, with convective storms dominating late spring-early summer months [54].The growing season lasts 180 days on average.The watershed is situated in the Southern Iowa Drift Plain, which has a rolling topography and many small streams.Hillslope gradients at the study sites range between 4% and 9%.Mollisols and Alfisols are the dominant soil orders found within the watershed.Loess deposits in the watershed, with depths up to 9 m [55], lead to very productive agricultural lands, but the soils are also highly erodible under the intensive land management [12,56].Three study sites under corn-soybean rotations with similar soil textures and surface gradients were carefully selected to isolate the role of tillage and row orientation on aggregate size distribution and stability (Figure 2; Table 1).Historically, these sites have been under agricultural production for approximately 100 years and were prairie lands before that [19].Although all the study sites were in a corn-soybean rotation, each had a unique management treatment and tillage intensity common to the watershed and region.Figure 3 provides a conceptual sketch of the three treatments.The first treatment was a 3-year, spring-till corn, ridge-till corn, and no-till soybean rotation.This treatment was a low-intensity tillage practice.The rows had an average ridge height of 10 cm (i.e., from crest to trough) and oriented along elevation contours, perpendicular to the downslope flow pathway (Figure 3a).This treatment is designated hereafter as Contour Ridge Tillage-Low Intensity (CRT-LI).The second treatment was similar to the first in terms of rotation and intensity; however, the rows were oriented parallel to the downslope flow pathway (Figure 3b).It is designated as Parallel Ridge Tillage-Low Intensity (PRT-LI).The third treatment was a spring-till corn, fall-till soybean rotation.The rows were oriented along elevation contours, perpendicular to the downslope flow pathway with a low roughness height, i.e., approximately 2 cm, due to the absence of ridges (Figure 3c).The tillage was more intense than CRT-LI and PRT-LI.This treatment is designated hereafter as Contour Tillage-High Intensity (CT-HI).

Flow Pathway Determination and Soil Sampling
The dominant flow pathways within each study field were determined by conducting geospatial analysis with the ArcMap 10.1 software (ESRI, Redlands, CA, USA) [62].LiDAR data obtained from Iowa Geodata [63] were used to develop a 1-m resolution raster digital elevation model (DEM) of each study area.The DEM was then filled using the spatial analyst toolset to remove any artifacts in the DEM.To calculate the direction that water would flow out of each cell, the spatial analyst toolset was used to determine the steepest descent path of each cell to one of eight neighboring adjacent cells based on elevation [64].From the steepest descents, a representative flow pathway was identified at each study site and selected for the sampling.As an example, the identified flow pathways for the CT-HI site are shown in Figure 4.
Four sampling locations were identified along each flowpath for each study site.At each The Soil Disturbance Rating (SDR) [57] was employed to characterize the level of tillage intensity for the different treatments.The SDR index considers the impact of the agricultural equipment used in a particular management strategy through soil inversion, mixing, lifting, shattering, aeration, and compaction treatments.CRT-LI and PRT-LI treatments had an SDR index value of 32, as opposed to a value of 53 for CT-HI [58].

Flow Pathway Determination and Soil Sampling
The dominant flow pathways within each study field were determined by conducting geospatial analysis with the ArcMap 10.1 software (ESRI, Redlands, CA, USA) [62].LiDAR data obtained from Iowa Geodata [63] were used to develop a 1-m resolution raster digital elevation model (DEM) of each study area.The DEM was then filled using the spatial analyst toolset to remove any artifacts in the DEM.To calculate the direction that water would flow out of each cell, the spatial analyst toolset was used to determine the steepest descent path of each cell to one of eight neighboring adjacent cells based on elevation [64].From the steepest descents, a representative flow pathway was identified at each study site and selected for the sampling.As an example, the identified flow pathways for the CT-HI site are shown in Figure 4. (see Figures 2 and 4).On collection, samples were maintained at field-moist conditions for later processing and analyses.The samples were collected in the spring time of 2014 and 2015 before any tillage disturbance or planting.

Aggregate Size Distribution
Figure 5A provides an overview of the sample processing steps.The collected soil samples were passed through an 8-mm sieve to remove coarse litter fractions, and then allowed to air dry for a period of three weeks [65].Approximately 125 g of the air-dried sample were dry sieved using an orbital shaker table and a nest of sieves with square openings.The following nine size classes were obtained: >2.00 mm; 2.00-1.Four sampling locations were identified along each flowpath for each study site.At each location, five replicate surface soil samples were collected from a depth of 0-5 cm within a 4-m 2 area (see Figures 2 and 4).On collection, samples were maintained at field-moist conditions for later processing and analyses.The samples were collected in the spring time of 2014 and 2015 before any tillage disturbance or planting.

Aggregate Size Distribution
Figure 5A provides an overview of the sample processing steps.The collected soil samples were passed through an 8-mm sieve to remove coarse litter fractions, and then allowed to air dry for a period of three weeks [65].Approximately 125 g of the air-dried sample were dry sieved using an orbital shaker table and a nest of sieves with square openings.The following nine size classes were obtained: >2.00 mm; 2.00-1.The weights of each size fraction were used to determine the sample dry mean weight diameter (DMWD) using the following equation [24,68]: The weights of each size fraction were used to determine the sample dry mean weight diameter (DMWD) using the following equation [24,68]: where x i is the mean diameter of each size class (mm); w i is the weight percentage of each size fraction with respect to the total sample; i is the respective size class; and n is the number of size classes.The dry mean weight diameter is an index that has been used to reflect degrees of tillage intensity on aggregate size distribution [24,69], and more importantly, the aggregate properties, i.e., more strongly aggregated soil will have a higher mean weight diameter, and vice versa.
After developing the dry aggregate size distribution, the size fractions within the SMAGG range (0.25-2.00 mm) were recombined and used for the aggregate stability tests.

Aggregate Stability Analysis
To determine aggregate stability under raindrop impact, air dried SMAGG samples were distributed evenly atop 0.25 mm sieves, which overlaid catch pans to collect slaked material during rainfall simulator tests (Figure 5B) [6].No pre-wetting of samples was performed to ensure the initial moisture content was similar for all samples [13, 70,71].A Norton Ladder Multiple Intensity Rainfall Simulator was used to mimic the effect of raindrop impact on SMAGG stability [72].The rainfall simulator was calibrated using a disdrometer to ensure a natural raindrop distribution was delivered over the sample area.The drop size distribution was controlled by the V-jet nozzle type and pressure setting (pressure was set at 6 PSI (pounds per square inch) and measured with the disdrometer [73].The raindrop size distribution agreed well with the Marshall-Palmer distribution, a commonly accepted distribution for natural raindrop sizes [74].
For each stability test, a rainfall intensity of 60 mm/h was used.This intensity corresponds to 1-h/25-year intensity storm [75] and it was selected based on prior studies that have examined the breakdown of soil aggregates under applied rainfall [15,32,76,77].To confirm that simulated rainfall was delivered uniformly, six testing locations (Figure 5B) were established underneath the simulator, which covered an effective area of 3 m 2 .Rainfall accumulation rates were quantified using the area of the collection containers (324 cm 2 ).The mean rainfall intensity for the six testing locations was 60.79 ± 1.38 mm h −1 , which corresponds to a variability of only 2% between locations.
Following established procedures [32], the rainfall was applied for 5 min for each test, supplying approximately 0.72 J of kinetic energy.The kinetic energy corresponding to this distribution of raindrops was estimated using the equations from [78][79][80].Following the 5-min test, both the material retained atop the 0.25 mm sieves and the material slaked into the catch pans were transferred to separate drying pans (Figure 5C).The retained and slaked fractions were then oven-dried at 60 • C for 24 h.The SMAGG stability, hereafter SMAGG STAB , was determined using the following relation [32]: where M r is the dry mass of the material retained on the sieve; M tot is the initial dry mass of the SMAGG sample; and M slk is the dry mass of material collected in the catch pan.

Quantification of Hydraulic Connectivity
Hydraulic connectivity was quantified by determining key overland flow variables with an enhanced version of the Water Erosion Prediction Project (WEPP) upland erosion model [39] for each of the flowpaths described in Section 2.2.Model inputs included a description of the downslope gradient/curvature, soil texture, management, and climate.The downslope gradient/curvature was obtained from the lidar-derived DEM described in Section 2.2.The soil texture along the flow pathway was specified from the sample analyses, as well as from previous studies at the sites [12,81] for deeper soil layers.The specified management practices corresponded to those described in Section 2.1.For consistency with the aggregate stability tests, the simulated storm event corresponded to the 25-year return period event (see Section 2.4).This simulated storm, which was obtained from rain gauge data [82], is representative of the study region and falls within the 90th-99th percentile of storms at the area [12,82].
Model simulation outputs included downslope profiles of the key overland flow variables that control the hydraulic forcing on surface aggregates [21], i.e., flow discharge, velocity, depth, and friction factor.These profiles were encapsulated into a single parameter profile, i.e., the applied shear stress profile, calculated and used as a proxy for hydraulic connectivity, for the sake of a physically-based representation of the hydraulic forcing on soil aggregates.The applied shear stress, τ (Pa), is related to the above key hydraulic parameters as follows: where f is the Darcy-Weisbach friction factor (-), u is the depth-averaged flow velocity (m/s) q is the unit flow discharge (m 3 /s/m), h is the flow depth (m), and ρ is the water density (kg m −3 ).To account for the soil resistance characteristics at the study sites and the resulting hydraulic forcing on dislodged, entrained, and transported aggregates by overland flow, a profile of the excess shear stress, δ, defined here as the difference between the applied shear stress and the critical soil shear stress, τ c , (i.e., δ = τ − τ c ) was derived for each flow pathway and used for the analyses.τ c is representative of the soil resistance and defined as the strength provided by interparticle forces of attraction or repulsion acting at the microscopic level, including electrostatic forces, van der Waals forces, hydration forces and biological forces [83,84].A calibrated baseline value of 3.5 Pa for τ c was adopted for the sites, following previous work [12,85,86].This baseline value was adjusted for tillage type with the WEPP relationship that accounts for roughness and consolidation effects [86].The adjusted critical shear stress values for the low-intensity tillage sites (CRT-LI and PRT-LI) and the high-intensity tillage site (CT-HI) were estimated as 5.46 Pa and 3.84 Pa, respectively.Finally, the event-averaged local excess shear stress profile was determined, representing the characteristic excess shear stress experienced by surface aggregates at any given point along the pathway.This was done to account for the hydraulic forces acting on in situ aggregates and against their interparticle cohesion forces.Values of the excess shear stress at the soil sampling locations were extracted from the profiles and used to represent the hydraulic characteristics corresponding to the sampled soils.

Statistical Analysis
Analysis of Variance (ANOVA) tests were made to determine if there were statistically significant differences between different hillslope positions for each management treatment.To meet the assumptions of the ANOVA test, the Shapiro-Wilk test was used to assess normality and the Brown-Forsythe Test was used for equal variance testing.In addition, Fischer-Snedecor (F) tests of coincidence of regression equations were performed to assess if aggregate size and stability versus downslope distance and local excess shear stress trendlines obtained for various management conditions were significantly different [87].Statistical analysis was performed using Sigma Plot 13.0.Significance level was set at 0.05.

Aggregate Size Distribution and Stability
At all sites, DMWD and SMAGG STAB decreased with downslope distance along the key flow pathways following a linear trend (Table 2; Figure 6).For each landscape position, the coefficient of variation for the measured values of DMWD and SMAGG STAB ranged between 2.2% and 26.3%, with the lower mean variability values reported for CRT-LI, and higher mean variability values for CT-HI (Table 2).F-tests for coincidence of trendlines showed that the trendline for each treatment, was significantly different from the trendlines of the other treatments (p < 0.05) for both DMWD and SMAGG STAB .Since management type and intensity did not vary in the downslope direction, the observed differences in the trends of DMWD and SMAGG STAB along the flow pathways and between the three types of managements should reflect the effects of hydraulic connectivity and tillage.Soil surface roughness produced by the various tillage practices introduces a level of overland flow control within a field, which affects hydraulic connectivity [49], and, thus, the hydraulic forcing on surface aggregates [21].This is particularly true with regards to the row orientation relative to the Soil surface roughness produced by the various tillage practices introduces a level of overland flow control within a field, which affects hydraulic connectivity [49], and, thus, the hydraulic forcing on surface aggregates [21].This is particularly true with regards to the row orientation relative to the direction of the main flow pathway in a field (oriented roughness), which has implications to the redistribution and downslope delivery of eroded material along the flow pathway [88].For CRT-LI and CT-HI, where the rows were oriented perpendicular to the flow pathway, a sequence of "steps and pools" has been observed along the downslope [85].This pattern reduced hydraulic connectivity by retarding or even impeding runoff in certain locations, thereby potentially limiting the mobilization and transport of aggregates as well as the hydraulic forcing exerted on them (i.e., reduced breakdown) [89].Conversely, for PRT-LI, where the oriented roughness configuration was parallel to the main flow pathways, higher flow concentration compared to CRT-LI and CT-HI treatments has been observed [19], which is expected to promote hydraulic connectivity.Row orientation may qualitatively explain the lower gradients of the DMWD and SMAGG STAB trend lines for CRT-LI and CT-HI compared to PRT-LI.

The Role of Hydraulic Connectivity
To quantitatively investigate the role of hydraulic connectivity on the trends presented in the previous section for the different management scenarios, the DMWD and SMAGG STAB values were plotted against their corresponding local excess shear stress (δ) values on the landscape, computed with the Papanicolaou et al. [39] framework (Figure 7).A consistent decrease was noted in both DMWD and SMAGG STAB with an increase in δ for all three management treatments.Two data groups were identified, i.e., low-intensity tillage (CRT-LI and PRT-LI), and high-intensity tillage (CT-HI), which fell into statistically different (p < 0.05) linear trend lines.Differences in the low-intensity and high-intensity tillage fields were apparent based on the gradient of the trendlines (Figure 7).The influence of tillage intensity on the observed trends is attributed to the profound effect that tillage operations can have on aggregate size through soil abrasion and soil quality through hydrologic weakening [6].The results indicate that a high SDR index value is associated with reduction in aggregate size.This likely leads to an increase in aggregate susceptibility, which results in decrease of the slope of both aggregate size and stability trendlines (Figure 7).This is in accordance with the findings from other studies that have shown that conventionally-tilled sites have smaller aggregate size fractions compared to reduced-or no-till sites [7,[90][91][92].Other studies have also observed higher aggregate stability in soils under reduced-tillage compared to conventional-tillage [26,93,94].The above suggest that both δ and SDR should be considered when examining the variability in DMWD and SMAGG STAB on a landscape in order to account for the roles of hydraulic connectivity and management.
Overall, the high Pearson correlation coefficients (Figure 7) suggest that the modeled local excess shear stress may be related to flowpath hydraulic connectivity and subsequent changes to aggregate size and stability for the different management conditions.The excess shear stress encapsulates the role of roughness features on aggregate size characteristics by accounting for the friction drag experienced by the aggregates.The estimated excess shear stress values exhibited different ranges over the flowpaths under the different management treatments, as illustrated in Figure 7.The excess shear stress for CRT-LI and CT-HI ranged from 0 to 1.10 Pa and 3.95 to 5.50 Pa, respectively, reflecting the runoff retardation effects of the "steps and pools".The lower excess shear stress values for CRT-LI compared to CT-HI were likely due to the larger ridge height of the CRT-LI, resulting in a larger retardation of flow and lower flow velocities.For the PRT-LI treatments, where accelerated runoff and sediment delivery have been observed due to flow concentration in the interrow areas [19], higher excess shear stress values were estimated compared to CRT-LI and CT-HI treatments, ranging from 0 to 8.16 Pa.These results support the notion that a lower degree of connectivity along a flow pathway arising from oriented roughness can reduce the variability in aggregate size and stability within a field as a result of lower absolute variability in the hydraulic forcing [19].Similar inferences may be drawn from other studies that have documented low variability of aggregate stability at different positions in pastures [95], where hydraulic connectivity is expected to be low, and, conversely, high variability of aggregate stability in steep and highly erodible landscapes [96], where hydraulic connectivity is expected to be high.

Relationship between Aggregate Size and Stability
The decreased SMAGGSTAB with increased local excess shear stress was remarkably similar to that of the DMWD (Figures 6 and 7).This observation prompted further investigation into the relation between DMWD and SMAGGSTAB, which revealed a positive relationship between the two aggregate characterization indices (Figure 8).F-tests for coincidence of regression lines showed that the trendlines for all sites were not statistically different (p < 0.05).The global regression model demonstrates that DMWD can explain 71% of the variability in SMAGGSTAB.
Although DMWD and SMAGGSTAB are important as standalone indices, the strong correlation between DMWD and SMAGGSTAB regardless of management suggests that the relationship between the two indices can be used as a simple predictor function for determining the aggregate stability from the aggregate size distribution at these sites and sites with similar attributes [19].This is advantageous when rainfall stability analyses are difficult to perform and only dry sieve analyses are possible.

Relationship between Aggregate Size and Stability
The decreased SMAGG STAB with increased local excess shear stress was remarkably similar to that of the DMWD (Figures 6 and 7).This observation prompted further investigation into the relation between DMWD and SMAGG STAB , which revealed a positive relationship between the two aggregate characterization indices (Figure 8).F-tests for coincidence of regression lines showed that the trendlines for all sites were not statistically different (p < 0.05).The global regression model demonstrates that DMWD can explain 71% of the variability in SMAGG STAB .
Although DMWD and SMAGG STAB are important as standalone indices, the strong correlation between DMWD and SMAGG STAB regardless of management suggests that the relationship between the two indices can be used as a simple predictor function for determining the aggregate stability from the aggregate size distribution at these sites and sites with similar attributes [19].This is advantageous when rainfall stability analyses are difficult to perform and only dry sieve analyses are possible.

Implications for Upland Erosion Modeling
The results of this study can be used to inform upland erosion models, by means of a better understanding and treatment of the physical nexus between management, surface hydrology, and aggregate breakdown, dislodgement, entrainment, and transport along the hillslope-scale continuum.Sidorchuk [97] conceptualized the soil detachment process as a conditional aggregate instability phenomenon, where detachment occurs when the driving hydrodynamic forces exceed the gravitational, hydrodynamic, geomechanical, interparticle and intraparticle resistance forces.A stochastic approach was considered by Sidorchuk [97], and it was demonstrated that the solution of the governing equations and the prediction of erosion rates in actual fields was intricate, mainly due to the lack of understanding and uncertainty in key soil properties such as aggregate size and stability.Even small variabilities in these properties were shown to introduce a high degree of uncertainty in erosion rates, rendering the theory implementation unfeasible at natural settings with our current knowledge and without the adoption of significant simplification assumptions.Results of this study give insight into the linkages between modeled or predicted hydraulic forcing and aggregates size and stability in the context of landscape connectivity.Specifically, they indicate an inverse relationship between the applied shear stress and the aggregate size and stability downslope along a flowpath.The results can guide further experimentation on the connection between these variables that could lead to the development of relationships that can be used in models to predict aggregate distribution and stability when the applied shear stress is known.Predicted aggregate stability and distribution in the models will govern material breakdown and mobilization, and subsequent transport by overland flow, leading to a better prediction of sediment, carbon and nutrient fluxes [19].
In addition, aggregate stability is a key parameter that controls the rainsplash and runoffinduced erosion across several spatial and temporal scales.De Roo et al. [98], and Laloy and Bielders [99] have successfully incorporated the aggregate stability in the mathematical formulation of erosion formulae.Legout et al. [87] developed a model of aggregate breakdown dynamics under rainfall in interrill areas using the mean weight diameter as an indicator of aggregate stability.Accurate

Implications for Upland Erosion Modeling
The results of this study can be used to inform upland erosion models, by means of a better understanding and treatment of the physical nexus between management, surface hydrology, and aggregate breakdown, dislodgement, entrainment, and transport along the hillslope-scale continuum.Sidorchuk [97] conceptualized the soil detachment process as a conditional aggregate instability phenomenon, where detachment occurs when the driving hydrodynamic forces exceed the gravitational, hydrodynamic, geomechanical, interparticle and intraparticle resistance forces.A stochastic approach was considered by Sidorchuk [97], and it was demonstrated that the solution of the governing equations and the prediction of erosion rates in actual fields was intricate, mainly due to the lack of understanding and uncertainty in key soil properties such as aggregate size and stability.Even small variabilities in these properties were shown to introduce a high degree of uncertainty in erosion rates, rendering the theory implementation unfeasible at natural settings with our current knowledge and without the adoption of significant simplification assumptions.Results of this study give insight into the linkages between modeled or predicted hydraulic forcing and aggregates size and stability in the context of landscape connectivity.Specifically, they indicate an inverse relationship between the applied shear stress and the aggregate size and stability downslope along a flowpath.The results can guide further experimentation on the connection between these variables that could lead to the development of relationships that can be used in models to predict aggregate distribution and stability when the applied shear stress is known.Predicted aggregate stability and distribution in the models will govern material breakdown and mobilization, and subsequent transport by overland flow, leading to a better prediction of sediment, carbon and nutrient fluxes [19].
In addition, aggregate stability is a key parameter that controls the rainsplash and runoff-induced erosion across several spatial and temporal scales.De Roo et al. [98], and Laloy and Bielders [99] have successfully incorporated the aggregate stability in the mathematical formulation of erosion formulae.Legout et al. [87] developed a model of aggregate breakdown dynamics under rainfall in interrill areas using the mean weight diameter as an indicator of aggregate stability.Accurate quantification of soil aggregate stability is therefore tied to the performance of erosion formulae and the accuracy of estimated sediment budgets.

Conclusions
The results of this study suggest that, along a flow pathway, management is not the only factor that affects aggregate proxies, but that hydraulic forcing is also a key factor that affects aggregate size and stability.Hence, hydraulic connectivity should be studied in conjunction with management for the evaluation of aggregate dynamics.Furthermore, within intensively managed systems, it may not be sufficient for DMWD and SMAGG STAB to be examined at soilscape scale, but rather at the landscape scale where the role of connectivity and transport can be adequately evaluated.
The low-intensity tillage sites examined in this study, CRT-LI and PRT-LI, had significantly higher DMWD and SMAGG STAB values than the high-intensity tillage site, CT-HI.At all sites, DMWD and SMAGG STAB values were found to decrease with increasing downslope distance and flow shear stress.The findings also suggested that the variabilities in the DMWD and SMAGG STAB values were related to the soil surface roughness that arose from management activities.Comparing contour and parallel low-intensity tillage sites, smaller variability in aggregate stability value was noted for the contour tillage site, reflecting the lower local excess shear stress experienced by surface aggregates stemming from lower hydraulic connectivity.Finally, a positive correspondence between SMAGG STAB and DMWD was observed for all the sites examined.A linear regression formula was developed with good performance, where DMWD explained 71% of the variability in SMAGG STAB .This relationship could be used as a simple function for predicting aggregate stability from the size distribution when aggregate stability measurements with rainfall simulators are difficult to perform.The relationship is empirical and site-specific, so it should be recalibrated for different settings, and more research is needed to establish a more universal relationship between aggregate size and stability.

Figure 2 .
Figure 2. Study sites in representative intensively managed agricultural fields of Clear Creek Watershed, IA.Management treatments include Contour Ridge Tillage-Low Intensity (CRT-LI), Parallel Ridge Tillage-Low Intensity (PRT-LI), and Contour Tillage-High Intensity (CT-HI).

Figure 2 .
Figure 2. Study sites in representative intensively managed agricultural fields of Clear Creek Watershed, IA.Management treatments include Contour Ridge Tillage-Low Intensity (CRT-LI), Parallel Ridge Tillage-Low Intensity (PRT-LI), and Contour Tillage-High Intensity (CT-HI).

Figure 3 .
Figure 3. Field pictures and conceptual sketches of the flow pathways sampled in the study sites: (a) Contour Ridge Tillage-Low Intensity (CRT-LI), a spring-till corn, ridge-till corn, no-till soybean lowintensity tillage system, where rows were oriented perpendicular to the downslope flow pathway (shown as red arrow); (b) Parallel Ridge Tillage-Low Intensity (PRT-LI), a spring-till corn, ridge-till corn, no-till soybean low-intensity tillage system, where rows were oriented parallel to the downslope flow pathway; and (c) Contour Tillage-High Intensity (CT-HI), a spring-till corn, fall-till soybean highintensity tillage system, where strip rows were oriented perpendicular to the downslope flow pathway.The yellow dashed lines in the field pictures represent the oriented roughness pattern.

Figure 3 .
Figure 3. Field pictures and conceptual sketches of the flow pathways sampled in the study sites: (a) Contour Ridge Tillage-Low Intensity (CRT-LI), a spring-till corn, ridge-till corn, no-till soybean low-intensity tillage system, where rows were oriented perpendicular to the downslope flow pathway (shown as red arrow); (b) Parallel Ridge Tillage-Low Intensity (PRT-LI), a spring-till corn, ridge-till corn, no-till soybean low-intensity tillage system, where rows were oriented parallel to the downslope flow pathway; and (c) Contour Tillage-High Intensity (CT-HI), a spring-till corn, fall-till soybean high-intensity tillage system, where strip rows were oriented perpendicular to the downslope flow pathway.The yellow dashed lines in the field pictures represent the oriented roughness pattern.

Figure 4 .
Figure 4. Sampling location selection using geospatial analysis: LIDAR data were used to generate a 1-m resolution digital elevation map for the study site.Elevation data were then analyzed to develop a flow direction map and isolate the dominant flow pathways (white arrow) that were sampled (green star) along the downslope.
Figure5Aprovides an overview of the sample processing steps.The collected soil samples were passed through an 8-mm sieve to remove coarse litter fractions, and then allowed to air dry for a period of three weeks[65].Approximately 125 g of the air-dried sample were dry sieved using an orbital shaker table and a nest of sieves with square openings.The following nine size classes were obtained: >2.00 mm; 2.00-1.70mm; 1.70-1.40mm; 1.40-1.18mm; 1.18-1.00mm; 1.00-0.85mm; 0.85-0.60mm; 0.60-0.25 mm; <0.25 mm.The sample was shaken at the lowest possible speed of 30 rpm and the mass of the size fraction retained on each sieve checked at 5 min and 10 mins to ensure that

Figure 4 .
Figure 4. Sampling location selection using geospatial analysis: LIDAR data were used to generate a 1-m resolution digital elevation map for the study site.Elevation data were then analyzed to develop a flow direction map and isolate the dominant flow pathways (white arrow) that were sampled (green star) along the downslope.

Figure 5 .
Figure5Aprovides an overview of the sample processing steps.The collected soil samples were passed through an 8-mm sieve to remove coarse litter fractions, and then allowed to air dry for a period of three weeks[65].Approximately 125 g of the air-dried sample were dry sieved using an orbital shaker table and a nest of sieves with square openings.The following nine size classes were obtained: >2.00 mm; 2.00-1.70mm; 1.70-1.40mm; 1.40-1.18mm; 1.18-1.00mm; 1.00-0.85mm; 0.85-0.60mm; 0.60-0.25 mm; <0.25 mm.The sample was shaken at the lowest possible speed of 30 rpm and the mass of the size fraction retained on each sieve checked at 5 min and 10 mins to ensure that it remained the same, and that the samples were not being significantly impacted by the slow circular shaking motion of the shaker[6,66,67].Geosciences 2018, 8, x FOR PEER REVIEW 8 of 19

Figure 5 .
Figure 5. (A) sample pre-processing steps to determine the dry mean weight diameter (DMWD) and isolate the small macroaggregate (SMAGG) size fraction to be tested for stability using the rainfall simulator; (B) plan view sketch of the rainfall simulator setup, front view photo of the rainfall simulator setup and plan view of the sieve setup for the aggregate stability tests; (C) steps followed to quantify SMAGG stability (SMAGG STAB ).

*
Landscape positon has a value of 1 to 4 to reflect location sample was collected (1 being top and 4 the bottom of the flowpath).

Figure 6 .
Figure 6.(a) relationship between dry mean weight diameter (DMWD) and downslope flowpath distance (x); (b) relationship between small macroaggregate stability (SMAGGSTAB) and downslope flowpath distance (x); the lines represent the linear trendlines obtained for each management condition.

Figure 6 .
Figure 6.(a) relationship between dry mean weight diameter (DMWD) and downslope flowpath distance (x); (b) relationship between small macroaggregate stability (SMAGG STAB ) and downslope flowpath distance (x); the lines represent the linear trendlines obtained for each management condition.

Geosciences 2018, 8 ,Figure 7 .
Figure 7. (a) relationship between DMWD and local excess shear stress (δ); (b) relationship between SMAGGSTAB and local excess shear stress (δ); light and dark green lines represent fitted lines obtained for low tillage intensities (Soil Disturbance Rating = 32) and high tillage intensities (Soil Disturbance Rating = 53), respectively.The Pearson correlation coefficient (r) is also shown for each tillage intensity dataset.

Figure 7 .
Figure 7. (a) relationship between DMWD and local excess shear stress (δ); (b) relationship between SMAGG STAB and local excess shear stress (δ); light and dark green lines represent fitted lines obtained for low tillage intensities (Soil Disturbance Rating = 32) and high tillage intensities (Soil Disturbance Rating = 53), respectively.The Pearson correlation coefficient (r) is also shown for each tillage intensity dataset.

Figure 8 .
Figure 8. Relationship between small macroaggregate stability (SMAGG STAB ) and dry mean weight diameter (DMWD) for all acquired samples.

Table 1 .
General site characteristics.

Table 1 .
General site characteristics.

Table 2 .
Mean and standard deviation values of the dry mean weight diameter (DMWD) and small macroaggregate stability (SMAGG STAB ) measured along the flowpaths for each sampled site.
* Landscape positon has a value of 1 to 4 to reflect location sample was collected (1 being top and 4 the bottom of the flowpath).

Table 2 .
Mean and standard deviation values of the dry mean weight diameter (DMWD) and small macroaggregate stability (SMAGGSTAB) measured along the flowpaths for each sampled site.