Microtopographic Controls on Erosion and Deposition of a Rilled Hillslope in Eastern Tennessee, USA

: Topography plays an important role in shaping the patterns of sediment erosion and deposition of different landscapes. Studies have investigated the role of topography at basin scales, whereas little work has been conducted on hillslopes, partially due to the lack of high-resolution topographic data. We monitored detailed topographic changes of a rilled hillslope in the southeastern United States using terrestrial laser scanning and investigated the inﬂuences of various microtopographic factors on erosion and deposition. The results suggest that the contributing area is the most important factor for both rill erosion and deposition. Rills with large contributing areas tend to have high erosion and deposition. Slope is positively related to erosion but negatively related to deposition. Roughness, on the other hand, is positively related to deposition but negatively related to erosion. Higher erosion and lower deposition likely occur on north-facing aspects, possibly because of higher soil moisture resulting from less received solar insolation. Similarly, soil moisture is likely higher in areas with higher terrain wetness index values, leading to higher erosion. This work provides important insight into the sediment dynamic and its microtopographic controls on hillslopes.


Introduction
Topography plays an important role in water-induced soil erosion and deposition, mainly through its impacts on surface flow pattern, velocity, and erosive power. Topography is also associated with soil properties, such as soil moisture, organic matter, particle size, affecting infiltration, runoff generation, and ponding [1]. Multiple topographic factors, such as slope length, gradient, and curvature, have been incorporated into various soil erosion models, including the Universal Soil Loss Equation (USLE) [2], Revised Universal Soil Loss Equation (RUSLE) [3], and Water Erosion Prediction Project (WEPP) [4]. However, topography has been primarily represented by a single factor in studies related to hillslope/experimental plots without the consideration of spatially varied terrain characteristics [5,6]. This single topography-related factor has been mainly used to emphasize slope length, gradient, or concavity/convexity [2][3][4][5], while topographic variation (e.g., roughness) has been often incorporated into other nontopographic parameters, for example, the cover factor in USLE.
It is generally accepted that topography affects the patterns of erosion and deposition. However, most of the emphasis has been focused on the basin spatial scale [7][8][9], with only a few attempts at plot or hillslope scales [10][11][12][13]. This may be due to a lack of the means to obtain high-resolution topographic data. Most publicly available digital elevation models (DEMs), particularly those obtained by airborne light detection and ranging (LiDAR), are at point densities or pixel resolutions of 50 cm to a few meters, which is not ideal for the

Study Area
Our study site is located in Loudon County, Tennessee (35 • 37 32.52 N, 84 • 12 59.69 W, Figure 1). This area is within the Ridge and Valley physiographic province on a terrace of the Little Tennessee River that originates from the Blue Ridge of the southern Appalachian Mountains and meanders westward into the Ridge and Valley. The climate is classified as the Köppen humid subtropical (Cfa) zone that is characterized by perennial precipitation and strong seasonality of hot summers and mild winters. High-intensity precipitation events are frequent due to tropical storms and mid-latitude cyclones. Based on U.S. climate data, the mean annual total precipitation is 1300 mm, and the mean annual temperature is 15 • C in this area (http://www.usclimatedata.com/climate/lenoir-city/tennessee/unitedstates/ustn0284, accessed on 5 March 2022). The study area is currently dominated by croplands and grassland, with chestnut oak (Quercus prinus) and eastern white pine (Pinus strobus) mixed forests on high elevation ridges and white oak (Quercus alba L.) and tulip poplar (Liriodendron tulipifera) associations in low-elevation valleys.
We monitored the topographic change of a 27 • engineered hillslope that was created during the construction of the Christenson Yacht facility in 2007 ( Figure 1). The soil of the hillslope is a gradually weathered silt and clayey Ultisols on shale parent material that is susceptible to rilling and gullying [23]. Particle size analyses of six soil samples collected from this site show consistency with the county soil survey with a soil texture of 70.0% clay, 14.7% silt, and 15.3% sand. These edaphic conditions are conducive to intense water erosion, and soil loss in this area has been impairing agricultural productivity and fragmenting the landscapes [24,25].
The hillslope faces southwest at an azimuth of 257° ( Figure 1). The altitude of the hillslope ranges from approximately 255 to 263 m above sea level. The hillslope has been almost without vegetation cover and has undergone intensive erosion. Grasses were later planted at the bottom and on the terraces above the hillslope for erosion control. Despite such efforts, rills with lengths up to 20 m have emerged on the hillslope surface from top to toe, detaching and transporting sediment that was subsequently deposited at the base of the hillslope. This study focuses on a vegetation-free section of the hillslope, extending about 18 m in length by 14 m in width. Because no manmade boundaries (earthen berms, metal sheets, etc.) were constructed, the hydrological process is not disturbed in this site, and the erodible materials are not exhaustive, making it ideal for long-term erosion monitoring and observation [26].

Data Acquisition
Data acquisitions for this study were previously explained in Lu et al. [10,11]. We conducted five field surveys on 10 December 2014, 7 March 2015, 11 June 2015, 13 September 2015, and 16 December 2015, respectively, using a 1550 nm wavelength FARO Fo- The hillslope faces southwest at an azimuth of 257 • (Figure 1). The altitude of the hillslope ranges from approximately 255 to 263 m above sea level. The hillslope has been almost without vegetation cover and has undergone intensive erosion. Grasses were later planted at the bottom and on the terraces above the hillslope for erosion control. Despite such efforts, rills with lengths up to 20 m have emerged on the hillslope surface from top to toe, detaching and transporting sediment that was subsequently deposited at the base of the hillslope. This study focuses on a vegetation-free section of the hillslope, extending about 18 m in length by 14 m in width. Because no manmade boundaries (earthen berms, metal sheets, etc.) were constructed, the hydrological process is not disturbed in this site, and the erodible materials are not exhaustive, making it ideal for long-term erosion monitoring and observation [26].

Data Acquisition
Data acquisitions for this study were previously explained in Lu et al. [10,11]. We conducted five field surveys on 10 December 2014, 7 March 2015, 11 June 2015, 13 September 2015, and 16 December 2015, respectively, using a 1550 nm wavelength FARO Focus3D X 330 TLS with an integrated global positioning system (GPS) and a dual axis compensator to keep data acquisition level and stable. We mounted the scanner on a tripod and set the height of the scanner at approximately 2 m above the ground. This TLS can scan the surrounding areas through a preset or customizable 360 • horizontal and 307 • vertical view at a nominal ranging radius of 330 m at 90% reflectance. The ranging accuracy of the scanner is ±2 mm at 50 m distance. During each scan, the scanner emits up to 976,000 laser pulses per second. The size of a laser spot is 2.25 mm at sensor exit and the laser beam diverges at 0.19 mrad (0.011 • ) as it moves further away. In our case, the ranging distance between the scanner and the hillslope surface is approximately 10 to 20 m, corresponding to a laser spot size of 4.15-6.05 mm [27]. We placed five identical spherical reference targets (ATS Scan Reference System with a diameter of 139 mm) around the hillslope before each survey (Figure 1). The targets are important for intra-survey registration, and we distributed them uniformly while avoiding any potential disturbance on the hillslope. We did not install permanent targets because of the concern of theft or vandalism. Because the scan targets are placed slightly differently for each survey, we only used the targets for registering intra-survey scans.
A Trimble GeoExplorer 6000 Series GeoXH™ handheld differential GPS (dGPS) was used to record the locations where the scanner and targets were placed. This dGPS unit has a built-in satellite-based augmentation system to perform the real-time correction on the coordinates once the system locks onto the strongest satellite signal. The nearest GPS base station is at the McGhee Tyson Air National Guard Base (35.81 • N, 84.00 • W), which is approximately 28 km away from our study site. We performed the dGPS measurement for at least 20 min at each location, allowing for the carrier to update the horizontal accuracy to 10 cm ± 2 parts per million (ppm) and the vertical accuracy to 20 cm ± 2 ppm. Since such level of accuracy is not sufficient for the cm level georeferencing, the coordinates collected by dGPS were only used to coarsely place the scans at their relative locations in the 3D virtual environment before registration. Further details of intra-survey and between-survey registrations are given in Section 3.2.
All TLS surveys were conducted on days with mild temperature and clear sky condition, and at least two days after any precipitation event. Each survey started with a coarse resolution (approximately 5 cm point spacing at 50 m radius) panoramic scan at the maximum viewing angle of the scanner unit (360 • horizontal by 307 • vertical) about 10 m away from the foot of the hillslope. Three scans were performed subsequently from three different locations at a similar distance from the hillslope. The subsequent scans were parameterized at a much finer resolution (approximately 1 cm spacing at a 50 m distance) to only scan the hillslope. For improved precision and comparable geometry of occlusion we tried to re-occupy the same scanning locations in each survey [28]. However, it was difficult to maintain identical locations due to the continuous change of the environment [29].

Preprocessing and DEM Generation
A relatively uniform density point cloud can be obtained by registering the point clouds produced by the TLS scans for each survey period (intra-survey registration). The registration was performed in the FARO SCENE ® software (version 5.1). This software provides a target-based registration that allows for intra-survey registration between different scans. An iterative closest point algorithm [30] is incorporated in the software to help align the scans on a target-to-target basis. We also visually checked the results of the registration and manually placed different scans using the scan targets and other fixed features, such as the wall of the yacht facility, in each scan to further improve the registration (Table 1).
A panoramic scan was conducted at a relatively coarser resolution to help the registration of the point clouds collected at different survey periods (henceforth called the between-survey registration). The first survey (14 December 2014) was used as the reference point cloud for the four between-survey registrations ( Table 1). The root-mean-square error (RMSE) of control points was used to assess the quality of the intra-survey and betweensurvey registrations ( Table 1). The RMSE is calculated in three-dimensional x, y, and z space. The RMSE of the intra-survey registration was calculated based on the offset of the five spherical targets, while the RMSE of the between-survey registration was calculated based on the offset of permanent artificial features (e.g., the walls of the yacht facility) captured in each survey.  [12,13,29,31].
We then rasterized the registered point cloud to produce a DEM for each survey period. No vegetation removal was necessary because this section of the hillslope was free of vegetation. After removing some outlier points, we converted each point cloud to a 1 cm pixel resolution DEM using a bilinear interpolation method in ArcGIS. The 1 cm pixel resolution was determined to best represent sediment estimates based on a previous study [10]. The bilinear interpolation algorithm determines the value of a grid cell using the linear combination of four neighboring cells [32]. It is suitable for continuous surface interpolation without distinct boundaries and has been widely used in converting LiDAR point clouds to raster DEMs [33][34][35].
All five DEMs were output to a local Cartesian projection. This projection focuses on the center point of the area of interest and defines a plane that is tangent to the spheroid at the center point, without accounting for the differences in z between corresponding points on the plane and the spheroid [36]. Note that the distortion in z is negligible due to the small extent of the study area. This projection best preserves the original quality of the data because the scanner stores the points first in polar coordinates that are then converted to Cartesian space.

DEM of Difference
We used the DEM of difference (DoD) method to calculate the erosion and deposition on the hillslope [37]. The DoD method has been widely used for geomorphic change detection where the differences between two or more DEMs obtained at different times can be quantified based on the following equation [37]: where t 0 is the initial time when the DEM was collected, t 1 is the subsequent time of data collection, and ε is the error. The calculation was performed on a cell-to-cell basis to produce a raster with cells of zero, positive, and negative values that spatially represent stable, depositional, and erosional areas on the hillslope, respectively. The error (ε) is related to the ranging error associated with the scanner (±2 mm), the registration errors (Table 1), the noise created by local topographic variation, and the error associated with environmental factors (atmospheric conditions, soil moisture, etc.). We used the Geometric Change Detection 6.0 (GCD) software with the capability to account for spatially varied uncertainties for the DoD calculation [37]. The GCD uses a fuzzy inference system that applies the fuzzy membership function to assign the uncertainties of each individual cell based on a set of input variables, including point density, roughness, and slope. For example, cells of high point density, low roughness, and low slope gradient are considered with low uncertainty and vice versa. When performing the DoD calculation, the elevation uncertainty is propagated along with uncertainties from various sources (e.g., registration, the survey method), and the elevation change at any cell is only included when the value is significantly (p < 0.05 or 95% confidence level) higher than the propagated uncertainty. A more detailed instruction on the GCD method is available in [37]. We repeated the DoD analysis across the time series, producing four DoD rasters from the five DEMs [10,11]. The elevation changes were only considered at the 95% confidence level. These four DoD rasters allow for the spatial mapping and estimation of both the volume and mass of erosion and deposition within this section of the hillslope [10,11].
Elevation is the altitude of a DEM cell, slope is the gradient, and aspect is the azimuth of the slope. Since aspect is a directional variable ranging from 0 • to 360 • , we transformed it to cosine (CosA) and sine (SinA) values for the analyses.
The roughness index (RI) describes topographic variation and irregularities in the elevation or standardized height of surface features for individual pixels. It is defined as where x 00 is the elevation of the central cell, and x ij is the elevation of each neighboring cell in a 3 × 3 window. A high roughness surface reduces erosion by increasing hydraulic resistance, dissipating the flow energy [39,40], and reducing runoff generation [41,42].
Contributing area (CA) is the upstream watershed area of a specific location on the hillslope. CA can be derived using flow direction, flow accumulation, and watershed analysis tools in ArcGIS. It has been widely used for channel network analysis [43,44] and soil erosion prediction [44].
Rill density (RD) is defined as the total length of rills per unit area [45]. It varies for slope steepness, slope length, runoff, soil, and precipitation conditions [46,47]. Channel depth (CD) measures the elevation difference between a specific cell location and its closest drainage divide (watershed boundary). Low CD values represent interfluve areas, while high values represent rill channels. This index has been used in soil classification and geomorphic feature analyses [48,49].
Topographic wetness index (TWI) was proposed to quantify the topographic control on hydrological processes [50]. TWI is defined as where CA is the contributing area, and θ is the slope angle in degrees. TWI is related to soil depth, infiltration rate, runoff generation, soil texture, and soil aggregates [51][52][53][54][55]. Slope length-gradient (LS) factor is a topographic index used in USLE and RUSLE for soil erosion prediction [2,3,7,56,57]. We derive the LS factor based on the method proposed in [7]: where A i,j−in is the contributing area at the inlet of grid cell (i,j); D is the cell size; and G x and G y represent gradient on x and y directions. The coefficient m is the length exponent of the USLE LS factor, representing the ratio of rill and interrill erosion. It is calculated by the β-value [57]: The β-value ranges between 0 (the ratio of rill to interrill erosion close to 0) and 1 and is defined as follows [57]: where θ is the slope angle. We derived the above nine topographic factors using ArcGIS 10 and SAGA GIS tools for the four 1 cm DEMs (December 2014, March 2015, June 2015, and September 2015). Each of the four DEMs was considered as the initial stages to which the DEM produced in the immediately subsequent survey was compared. The topographic changes occurred between each survey and the next are produced using the DoD method described in Section 3.3.
We used the "sub-watershed" method described in [11] to derive summary statistics (median, sum, or percentage value) for each factor, as well as the erosion and deposition amounts along the rills. The rill network was derived using the hydrology tools in Tau-DEM [58]. A rill basin refers to all upstream pixels contributing to the lowest outlet point along a rill channel. We generated a set of outlet points at a 0.05 m interval along the rills to segment a rill basin into a series of non-overlapping polygons (sub-watersheds). Although the generated sub-watersheds have various shapes and sizes, a fixed interval guarantees a uniform sampling distance along the rill basin. Please check the reference [11] for detailed description of the "sub-watershed" method.

Correlation and Regression Analysis
Pearson's correlation analysis was used to examine the relationships among erosion, deposition, and topographic factors. This analysis produces a correlation coefficient (r) to measure the strength and a p-value to represent the statistical significance of the correlation between each pair of two variables.
We used the linear regression models to examine the relationship between erosion/ deposition and topographic factors. The linear models were chosen because they provide clarity and interpretability compared to machine learning models, although they may compromise some predictability. A linear regression model can be expressed as follows: where y is the dependent variable; x 1 , x 2 , . . . , x n are the independent variables; β 0 , β 1 , . . . , β n are the regression coefficients; and ε is the error. Both ordinary least squares (OLS) and quantile regression (QR) were used in the analysis. OLS estimates the conditional mean of the dependent variable based on the estimation of regression coefficients (β) using the minimum of the sum of squared residuals: OLS is usually used when both dependent and independent variables are normally distributed, and it is not very robust against datasets that have potential outliers.
QR is a modification of OLS that estimates the regression coefficient using the conditional quantile [60]:β where τ is the τth quantile (0 < τ < 1) that is defined as QR is more robust against the datasets with potential outliers and provides better performance when variables have a non-Gaussian distribution. QR is also ideal when the distribution of the data is skewed, especially when the variables do not exhibit homoscedasticity (similar variance for all independent variables) [60].
Spatial autocorrelation may exist among these topographic factors, violating the assumption of independence for common statistical practices [61]. In our case, various sources might contribute to the spatial autocorrelation of the topographic factors, for example, through hydrological and geomorphic processes associated with erosion and deposition on the hillslope. We used a group k-fold cross-validation method to minimize the influence of spatial autocorrelation while achieving the optimal quantile τ without overfitting. The group k-fold cross-validation method divides data into k groups based on a categorical variable; the cross-validation is repeated for k times, and for each iteration, k-1 groups are used to train the model and the one group that is left out is used to validate the model. The k results from the iterations are eventually averaged to yield an estimation. In this study, a set of parallel rills that were developing on the hillslope comprised our cross-validation dataset. To reduce spatial autocorrelation, we divided the entire dataset into different rill basins based on the rill drainage points close to the foot of the hillslope. We then ensured that the segments within the same rill basin cannot be included in both the training and testing sets.
Since we included various topographic factors that may or may not share redundant information, it is possible for our model to produce biased estimation unless variables are selected in an objective fashion. We used the recursive feature elimination (RFE) with external validation to select the fewest possible variables that preserve the best predictive ability [62]. Sometimes, a backward feature selection approach is used for statistical models: the model loads all variables at the beginning and the least important variable is dropped at each iteration until the model performance converges at an optimal set of variables. However, such an approach is susceptible to overfitting of the training dataset [63]. The RFE approach accounts for such issues by enclosing the backward feature selection within a resampling procedure. For each iteration, the data are partitioned into different training and testing subsets, and the performance of the model is evaluated using the mean squared error (MSE) of the testing subset. The variables included in the final model are determined based on a consensus ranking from all iterations [62].    Figure 2 illustrates the spatial pattern of erosion and deposition on this hillslope section during the four periods. Generally, the rills are more dynamic compared to interrill areas, with the majority of the detectable changes occurring either within or in proximity to the rills. The majority of the interrill areas are relatively stable with only a few detectable "hotspots" of erosion or deposition.   Figure 2 illustrates the spatial pattern of erosion and deposition on this hillslope section during the four periods. Generally, the rills are more dynamic compared to interrill areas, with the majority of the detectable changes occurring either within or in proximity to the rills. The majority of the interrill areas are relatively stable with only a few detectable "hotspots" of erosion or deposition.   Table 3 shows the results of Pearson's correlations between the dependent variables and independent variables. Both erosion and deposition are correlated with several in-dependent variables, and many independent variables are correlated with one another. High correlations are observed between some of the variables. For example, the correlation coefficient between slope and the LS factor is 0.83, and the correlation between slope and RI is 0.97. The LS factor is highly correlated with CA (0.73) and RI (0.84). The correlations among independent variables might lead to biased estimation of coefficients and affect the efficacy of the linear regression models. We used RFE to select the fewest variables that could achieve the minimal MSE. Table 4 shows the results of the RFE. The erosion model and the deposition model achieved the least MSEs with eight and nine variables, respectively. The variables in Table 4 are ranked in correspondence to their importance to the model. In the erosion model, CA is the most important variable. The model with only CA as the independent variable has a relatively high MSE of 0.286. With more variables included in the model, the MSE gradually decreases until it increases again. The second most important factor is slope, followed by CosA, CD, RI, TWI, LS, and elevation. For the deposition model, CA is also of the highest importance, followed by elevation, CD, CosA, TWI, LS, slope, RD, and RI.  Figure 3 shows the comparison between the observed and predicted erosion/deposition using OLS and QR models. The OLS models produced R 2 values of 0.19 and 0.31 for erosion and deposition, respectively. The QR models produce R 2 values of 0.29 and 0.34, respectively. The QR models have stronger explanatory capability compared to the OLS models, possibly because the amounts of erosion and deposition do not follow a normal distribution on the hillslope. The QR models also produced less-biased prediction at the extreme values when the optimal quantile parameter was adopted. The predictions made by the QR models showed a better agreement with the observed values compared to those of the OLS models (the black and red dashed lines are closer in QR models).  Tables 5-8 list the coefficients of RFE-selected topographic variables in the QR models for erosion and deposition during the four scanned periods, respectively, with the quantile (τ) that produces the smallest cross-validation MSE. Generally, CA is associated with higher erosion and deposition, although its coefficient (absolute value) is relatively higher in the erosion model, except for the period from June to September. A steeper slope is associated with higher erosion and lower deposition. In our study area, the hillslope has been incised with well-developed rills, and rill sidewalls are steeper than the interrill and rill floor areas. During rainfall events, erosion would be more significant on rill sidewalls than interrill and rill floor areas. The coefficient of slope, representing the impact of slope on the magnitude of erosion/deposition, is the highest during the period of March to June. This is likely because the soil surface is relatively loose at the beginning of this Tables 5-8 list the coefficients of RFE-selected topographic variables in the QR models for erosion and deposition during the four scanned periods, respectively, with the quantile (τ) that produces the smallest cross-validation MSE. Generally, CA is associated with higher erosion and deposition, although its coefficient (absolute value) is relatively higher in the erosion model, except for the period from June to September. A steeper slope is associated with higher erosion and lower deposition. In our study area, the hillslope has been incised with well-developed rills, and rill sidewalls are steeper than the interrill and rill floor areas. During rainfall events, erosion would be more significant on rill sidewalls than interrill and rill floor areas. The coefficient of slope, representing the impact of slope on the magnitude of erosion/deposition, is the highest during the period of March to June. This is likely because the soil surface is relatively loose at the beginning of this period after the winter freeze-thaw cycles so that the hillslope is more vulnerable to erosion. Similar to slope, a higher roughness is associated with higher deposition and lower erosion because rugged terrain surface tends to dissipate flow power, reducing transport capacity. Higher elevation leads to lower erosion and lower deposition, which is consistent with the field observation that the areas closer to the top of the hillslope are generally more inactive compared to the foot area. Higher CD values lead to less erosion and more deposition. On this hillslope, areas of high CD values are usually where rills are most deeply incised, such as rill floors; the deeper the rills grow, the more compact and less erodible materials (close to the bedrock) are.

Correlation and Regression Analyses
Several other topographic factors affect erosion and deposition through their associations with soil moisture content. Our results suggest that a larger CosA leads to higher erosion and lower deposition. CosA represents the aspect facing a north-south direction, and a larger CosA is closer to north. North-facing areas are more likely have higher soil moisture content in the Northern Hemisphere. These areas are easily saturated, leading to the generation of more runoff during a rainfall event. Similarly, the areas with higher TWI values have higher erosion and lower deposition because TWI estimates the likelihood of water accumulation in an area with elevation differences [51].

Seasonal Variation of Erosion and Deposition
The patterns of erosion and deposition are related to the different precipitation and temperature conditions during the four scan periods. Less erosion and more deposition occurred from December 2014 to March 2015 compared to other periods. This is likely caused by the relatively limited transport capacity during this period. Most days of this period had temperatures below 0 • C at this site, and the observed topographic change might be driven by the freeze-thaw process. Figure 4 illustrates the changes in daily high temperature, low temperature, precipitation, and snow at Lenoir City, TN, USA (35.79 • N, 84.26 • W, www.usclimatedata.com, accessed on 5 March 2022), approximately 18 km from the study site. The daily low temperature dropped below 0 • C starting late October and fluctuated above and below 0 • C for most days of November, December, and January. Meanwhile, precipitation events may keep the soil moist, allowing for freeze-thaw cycles with the diurnal variation of the temperature. Barnes et al. [64] found that a thin soil layer on the surface can be heaved by ice crystals formed during freezing seasons, and these crystals lead to failures and within-channel deposition during thawing. The mass failures along the sidewalls created by freeze-thaw cycles may widen cross-sections, reducing flow velocities within the rills [65]. The DoD result during the first period suggests that limited transport capacity created more in-channel deposition in the first period, leading to excessive sediment supply in the subsequent period (March to June 2015), shown as the high erosion rate.
As listed in Tables 5-8 and shown in Figure 5, the coefficients of the same topographic factor are different in the four survey periods. Note that negative coefficients represent positive influences, while positive coefficients represent negative influences for the topographic factors in the erosion models because erosion values are represented by negative elevation changes. The erosion models show that CA, slope, CosA, and TWI have the highest positive influences, whereas CD, LS, and elevation show the highest negative influences during the March-June period in these four periods. The highest positive or negative influences of these factors during the March-June period are likely associated with the excessive sediment supply following the winter period. Because of the excessive sediment supply, all factors related to surface flow (CA, LS, elevation, and slope), soil moisture (TWI and CosA), and erodibility (CD) actively act on sediment erosion and transportation during the March-June period. In contrast, the decreasing influences of these factors in the following periods likely reflect the limited sediment supply on hillslope after eroding and transporting the excessive sediments during the March-June period. In addition to the spring period (March-June 2015), our results also indicate that slope has a higher influence during the winter (December 2014-March 2015) period than the other two periods. This might suggest that slope, at the microtopographic scale, is strongly associated with erosion during freeze-thaw cycles as well.
For the deposition models, CA has the greatest influence during the September-December 2015 period and the least influence during the December 2014-March 2015 period. CosA and TWI are both negatively related to deposition in all deposition models, with the highest negative influences in the March-June 2015 period. These two factors are associated with soil moisture, and sediments are more erodible or transportable in higher soil moisture areas when sediment supply is excessive, consistent with the results that these two factors have the highest impacts on the erosion models during the same period. Slope had higher negative impacts during March-June 2015 (−0.227) and June-September 2015 (−0.210) periods, suggesting that the areas of higher slope gradients (usually rill headwalls and sidewalls) have less chance of deposition. RI, associated with the flow friction imposed by surface variation, has the highest influence on deposition during the winter period (December 2014-March 2015). This is likely caused by the freeze-thaw cycles, which create more surface irregularity and friction to surface flow, leading to more sediment deposition.

Factors Controlling Erosion and Deposition
Microtopographic factors show varied importance regarding their impact on erosion and deposition at different scales. CA is the most important variable in all models; both high erosion and deposition occur on the areas with high CA values, reflecting the influence of relative location on the hillslope and within a rill basin. High erosion and low deposition occur on the areas with larger TWI values and on north-facing aspects (larger CosA). Both TWI and aspect are related to soil moisture content, which is likely lower on the areas with lower TWI or CosA (south-facing) values. During a rainfall event, relatively dry soil can absorb more rainfall, postponing and reducing runoff generation. CD is another important factor in most models. CD represents the relative location along the rill cross-sectional profile. A larger CD value suggests low erosion and high deposition in the models. This is because rills have been well developed on the hillslope, and rill floors (high CD areas) have reached a compact layer where the parent materials are not fully decomposed. This compact layer is less erodible than the topsoil of the hillslope. On the other hand, sediments eroded from steep rill sidewalls are likely deposited, at least temporarily, on rill floors where CD is relatively high.

Limitations of This Work
The complex nature of the erosion and deposition processes introduces challenges to quantifying sediment dynamics at a hillslope scale using grid-based approaches [31,66,67]. The choice of surveying hardware, software, and data processing and analytical methods should be in correspondence with the size of the study area, the features to be monitored, the acceptable level of error, and other factors, such as the frequency of data acquisition and logistic availabilities [10]. TLS is useful in our study by allowing for rapid survey of the hillslope surface without direct interference of the erosional and depositional processes, and the obtained high-resolution topographic data enable the differentiation and quantification of the spatially varied erosion and deposition. However, limited by the uncertainty from various sources, it is still difficult to detect the microtopographic changes on interrill areas.
A closer field examination shows that the redistribution of sediment in our study area is further complicated by other factors, including surface armoring and crusting ( Figure 6). The erosive power gradually exhausts relatively finer particles (e.g., clay, silt, and fine sand), leaving behind coarser materials, such as some pebbles (Figure 6a). Although our model considered factors such as roughness and TWI as proxies of soil properties, including particle size, organic matter, and soil moisture, these factors are only approximations of the ground truth at best. In situ measurements are impossible as mass sampling at such fine resolution can be expensive, time-consuming, and would inevitably disturb the sediment on the hillslope. As listed in Tables 5-8 and shown in Figure 5, the coefficients of the same topographic factor are different in the four survey periods. Note that negative coefficients represent positive influences, while positive coefficients represent negative influences for the topographic factors in the erosion models because erosion values are represented by negative elevation changes. The erosion models show that CA, slope, CosA, and TWI have the highest positive influences, whereas CD, LS, and elevation show the highest negative influences during the March-June period in these four periods. The highest positive or negative influences of these factors during the March-June period are likely associated with the excessive sediment supply following the winter period. Because of the excessive sediment supply, all factors related to surface flow (CA, LS, elevation, and slope), soil moisture (TWI and CosA), and erodibility (CD) actively act on sediment erosion and transportation during the March-June period. In contrast, the decreasing influences of these factors in the following periods likely reflect the limited sediment supply on hillslope after eroding and transporting the excessive sediments during the March-June period. In addition to the spring period (March-June 2015), our results also indicate that slope has a higher influence during the winter (December 2014-March 2015) period than the other two periods. This might suggest that slope, at the microtopographic scale, is strongly associated with erosion during freeze-thaw cycles as well. CosA and TWI are both negatively related to deposition in all deposition models, with the highest negative influences in the March-June 2015 period. These two factors are associated with soil moisture, and sediments are more erodible or transportable in higher soil moisture areas when sediment supply is excessive, consistent with the results that these two factors have the highest impacts on the erosion models during the same period. Slope had higher negative impacts during March-June 2015 (−0.227) and June-September 2015 (−0.210) periods, suggesting that the areas of higher slope gradients (usually rill headwalls osition in the rills. We also observed surface crusting on the hillslope (Figure 6b), which can affect raindrop impact on interrill areas [68]. The sporadic crusting may create variation in erodibility between different interrill areas. It is also possible that microorganisms (worms) or animals may disturb the site between surveys, causing potential errors. In addition, the erosion and deposition processes are also affected by rainfall characteristics, such as intensity and duration. More work is needed in the future to examine the impacts of rainfall characteristics on hillslope erosion and deposition. Other factors, such as the pebbles commonly observed on the rill floor, also likely affect the erosion/deposition patterns by altering surface hydrology. The pebbles help prevent raindrop splash and sheet erosion on interrill areas by reducing raindrop impact and dissipating flow energy; they also protect rill floors from further entrenchment. The increase in roughness because of the pebbles reduces the transport capacity, leading to deposition in the rills. We also observed surface crusting on the hillslope (Figure 6b), which can affect raindrop impact on interrill areas [68]. The sporadic crusting may create variation in erodibility between different interrill areas. It is also possible that microorganisms (worms) or animals may disturb the site between surveys, causing potential errors. In addition, the erosion and deposition processes are also affected by rainfall characteristics, such as intensity and duration. More work is needed in the future to examine the impacts of rainfall characteristics on hillslope erosion and deposition.

Conclusions
This study examined the influences of microtopographic factors on sediment dynamics on a rilled hillslope in eastern Tennessee, USA, based on linear regression models. Generally, the QR model showed a stronger predictive ability compared to the OLS model, with 29% of the variability for erosion and 34% for deposition explained by the microtopographic variables. The coefficients of nine factors (contributing area, roughness, slope length, gradient, rill density, aspect, TWI, LS factor, and channel depth) reflect the level of their influences on the prediction of erosion or deposition in the QR models. A larger CA leads to higher erosion and deposition. A steeper slope increases erosion and reduces deposition. Roughness is positively related to deposition and negatively related to erosion. Rill floors with larger depth values tend to have less erosion and more deposition. The cosine component of the aspect is associated with higher erosion and lower deposition, possibly due to higher soil moisture on the north-facing aspects with lower solar insolation received. TWI, another index positively associated with soil moisture, is positively related to erosion, and negatively related to deposition. The influences of each topographic factor on erosion and deposition vary in the four periods. The highest positive or negative influences of most topographic factors on erosion during spring are likely associated with the excessive sediment supply accumulated in winter. Then, after eroding and transporting the excessive sediments in spring, the influences of these factors decrease in the following summer and fall periods. The influences of these topographic factors on deposition are mostly consistent with (opposite to) the influences on erosion.

Data Availability Statement:
The data related to this study are available on request from the authors.