Seasonal Precipitation Variability and Gully Erosion in Southeastern USA

: This study examines the relationship between gully erosion in channels, sidewalls, and interﬂuves, and precipitation parameters (duration, total accumulation, average intensity, and maximum intensity) annually and seasonally to determine seasonal drivers for precipitation-related erosion. Ordinary Least Square regression models of erosion using precipitation and antecedent precipitation at weekly lags of up to twelve weeks were developed for three erosion variables for each of three geomorphic areas: channels, interﬂuves, and sidewalls (nine models in total). Erosion was most pronounced in winter months, followed by spring, indicating the inﬂuence of high-intensity precipitation from frontal systems and repeated freeze-thaw cycles in winter; erosion in summer was driven by high-intensity precipitation from convectional storms. Annually, duration was the most important driver for erosion, however, during winter and summer months, precipitation intensity was dominant. Seasonal models retained average and maximum precipitation as drivers for erosion in winter months (dominated by frontal systems), and retained maximum precipitation intensity as a driver for erosion in summer months (dominated by convectional storms). In channels, precipitation duration was the dominant driver for erosion due to runo ﬀ -related erosion, while in sidewalls and interﬂuves intensity parameters were equally important as duration, likely related to rain splash erosion. These results show that the character of precipitation, which varies seasonally, is an important driver for gully erosion and that studies of precipitation-driven erosion should consider partitioning data by season to identify these drivers. with coe ﬃ cients of determination ranging from R 2 = 0.245 to R 2 = 0.49 (except for IErosion in summer at R 2 = 0.131 and SDep in winter at R 2 = 0.087), suggesting that precipitation is an important driver for erosion in these months, no matter the metric used. Moreover, these results show that the character of the precipitation is an important driver for erosion; antecedent precipitation has an inﬂuence on erosion in the following weeks and months and it varies with season. -related erosion, while in sidewalls and interﬂuves, intensity parameters were equally important as duration, likely related to rain splash erosion. This research shows that soil erosion is seasonally variable and an understanding of the seasonal pattern of soil erosion with respect to precipitation-related drivers improves the potential to achieve strategic conservation measures.


Introduction
Gully erosion is a global problem, particularly in the southeastern United States, where erodible soils, high relief, and climatic and meteorological factors encourage soil erosion. Gully erosion is one of the most dangerous forms of soil degradation, which is caused by natural and anthropogenic activities. Gullies are composed of several continuous or discontinuous channels and rills with varying slopes, which may later develop into deep trenches, inhibiting effective remediation by tillage. Gully erosion can initiate from anthropogenic factors like farming or grazing on susceptible soils, increased runoff from land-use changes due to logging or construction, and poor vegetative cover from wildfire or high soil salinity. Additionally, natural drivers for soil erosion are meteorological variables, topography, and soil type and texture [1,2].
Changes in land use can increase soil erosion. Vast regions of the United States experienced soil erosion when forested lands were converted to croplands in the late 19th century and the early 20th century [3]. Estimates of the volume of soil erosion in the United States caused by both sheet and rill erosion combined is 6.7 Mg/ha/y in cultivated cropland, 0.90 Mg/ha/y on federal lands, and 1.55 Mg/ha/y in pasture lands [3]. Considerable land area in the southeastern US was converted from forest to agriculture to support cotton farming in the 1800s and pasture for animal grazing [4,5].
Ridge physiographic province and consisted of northeast-southwest trending parallel limestone valleys (Maynardville Formation) and sandstone or shale ridges (Nolichucky Formation) [17]. The 1.5 ha study area was located on a grass and shrub hillslope surrounded by forest (on the ridges) and pasture (in the valleys). Soils were highly erodible fine-grained silt and clay Ultisols (Collegedale-Etowah complex (CeD3)) with an average erodibility factor (RUSLE K-factor) of 0.28, indicating susceptibility to raindrop impact and transport by surface runoff [18]. The region has a humid subtropical climate (Köppen Cfa) with year-round precipitation of 1070 mm (42 in) annually and an average annual temperature range from 1.1 °C (34 °F) in January to 23.3 °C (74 °F) in July. The National Oceanic and Atmospheric Administration describe Tennessee's winter precipitation as dominated by the polar front and summer precipitation that results from convectional systems. September and October are the driest months. A detailed description of the site setup can be found in [13,19] and is summarized as following. Steel erosion pins were installed in transects throughout the 100 m × 100 m gullied zone. Each transect spanned interfluves, sidewalls, and the gully channel to assess erosion in these three morphological settings. In total, 105 erosion pins were installed, 34 (1 m × 5 mm) pins in channels, and the remaining (0.5 m × 5 mm) pins in interfluves (29 pins) and sidewalls (42 pins). From 23 May 2012 to 22 August 2018, pin length was recorded approximately weekly for each pin using a folding ruler. Pin attrition occurred periodically over the study period, such that some pins were eroded, damaged, or dislodged by animals. Therefore, in May 2015, 43 new pins were installed and 3 damaged pins were replaced, bringing the total number of pins to 105. The nature of the site surface limited access during and immediately after rain events, and over the six-year period, pin length was recorded 294 times. The difference between the exposed lengths of each pin was calculated between one measurement period and the next, and this dataset of pin change was compared to precipitation data to identify important drivers for erosion in each morphological setting.
For each setting, we created three erosion variables: (1) average of the absolute value of change (Avg|Ch|); (2) average of only positive changes in pin lengths (deposition) from one measurement period to the next (AvgDep), and; (3) average of only negative changes in pin lengths (erosion) from one measurement period to the next (AvgErosion). In prior research, a fourth variable, average change, was generated, however, because of a balance of erosion and deposition, especially in channels, the average change remained near zero and was not a useful parameter to capture weekly A detailed description of the site setup can be found in [13,19] and is summarized as following. Steel erosion pins were installed in transects throughout the 100 m × 100 m gullied zone. Each transect spanned interfluves, sidewalls, and the gully channel to assess erosion in these three morphological settings. In total, 105 erosion pins were installed, 34 (1 m × 5 mm) pins in channels, and the remaining (0.5 m × 5 mm) pins in interfluves (29 pins) and sidewalls (42 pins). From 23 May 2012 to 22 August 2018, pin length was recorded approximately weekly for each pin using a folding ruler. Pin attrition occurred periodically over the study period, such that some pins were eroded, damaged, or dislodged by animals. Therefore, in May 2015, 43 new pins were installed and 3 damaged pins were replaced, bringing the total number of pins to 105. The nature of the site surface limited access during and immediately after rain events, and over the six-year period, pin length was recorded 294 times. The difference between the exposed lengths of each pin was calculated between one measurement period and the next, and this dataset of pin change was compared to precipitation data to identify important drivers for erosion in each morphological setting.
For each setting, we created three erosion variables: (1) average of the absolute value of change (Avg|Ch|); (2) average of only positive changes in pin lengths (deposition) from one measurement period to the next (AvgDep), and; (3) average of only negative changes in pin lengths (erosion) from one measurement period to the next (AvgErosion). In prior research, a fourth variable, average change, was generated, however, because of a balance of erosion and deposition, especially in channels, the average change remained near zero and was not a useful parameter to capture weekly and longer-term erosion on-site [13,14,[19][20][21]. Therefore, in this study, we have retained the three variables described above.
A Davis Vantage Pro wireless weather station (KTNJONES12, data available at https://www. wunderground.com/dashboard/pws/KTNJONES12) was located 350 m from the research site, and recorded precipitation, pressure, temperature, and wind data at five-minute intervals. Occasional data gaps were filled with data from a neighboring station 1.6 km away (KTNJONES7, data available at https://www.wunderground.com/dashboard/pws/KTNJONES7), with only 21 of 2282 study days missing weather data. See [19] for a detailed list of weather data gaps and coverage.
From these data, four precipitation parameters were generated for each measurement period: (1) Duration (total minutes of rainfall); (2) Total Accumulation (total precipitation in mm); (3) Average Intensity in mm/min (Total Accumulation/Duration), and; (4) Maximum Intensity in mm/min (the greatest station-reported rain rate during the measurement period). The rain rate is a smoothed function of rain accumulation over time that is calculated using the ratio of the tipping bucket depth-adjusted volume to the time between tips. As rainfall tapers off, the rate drops but does not reach zero immediately upon cessation of precipitation. Instead, it smooths the rate to more accurately represent how precipitation naturally tapers over an area at the end of a rain storm [22].
Prior research has shown that antecedent precipitation may be an important factor in erosion, and therefore a series of antecedent precipitation parameters were generated for the prior eleven measurement periods, for each of Duration, Total Accumulation (TotAcc), Average Intensity (AvgInt), and Maximum Intensity (MaxInt). These antecedent lagged variables were named Duration-1, Duration-2 ... Duration-11, TotAcc-1, TotAcc-2 . . . and so-on, a total of 48 precipitation parameters, which we refer to as lagged precipitation parameters.
The relationship between erosion variables and all precipitation parameters was assessed with Spearman correlation coefficients. Ordinary Least Squares (OLS) regression models were created for the nine erosion variables using the set of current and lagged precipitation parameters. Further, because seasonal variability in erosion was observed in prior studies [13,19], the data were partitioned by season: winter (December, January, February); spring (March, April, May); summer (June, July, August); and autumn (September, October, November). OLS regression models were generated for the erosion variables using the precipitation parameters for each of the seasonal datasets.

Precipitation
Precipitation accumulation for each measurement period had an annual mean of 22.2 mm, with the highest seasonal mean accumulation in winter (26.3 mm) and spring (24.3 mm), and the lowest in autumn (15.3 mm) (Table 1). Likewise, the duration of precipitation had an annual mean of 278.7 min, but the longest seasonal mean duration was received in winter (424.7 min), and the shortest in autumn (192.5 min). Both average and maximum precipitation intensity were higher in summer months (0.1 mm/min and 108.1 mm/min, respectively) compared to the annual values of these parameters (0.08 and 71.5 mm/min, respectively).
The study area experienced year-round precipitation, however, most of the accumulation was in winter (frontal systems) and summer (convectional storms) ( Figure 2). September and October were the driest months, most notably in 2012, 2013, and 2016. The most intense rains occurred in summer months, for example, see high values for Average Intensity (AvgInt) in the summer of 2012, 2014, 2016, 2017, 2018, and to a lesser degree 2013 and 2015. One may also notice that when Total Accumulation (TotAccum) was high and Duration was low, Maximum Intensity (MaxInt) was also high because it follows that higher intensity rainfall occurred when high rainfall totals were received in a short time-period. The study area experienced year-round precipitation, however, most of the accumulation was in winter (frontal systems) and summer (convectional storms) ( Figure 2). September and October were the driest months, most notably in 2012, 2013, and 2016. The most intense rains occurred in summer months, for example, see high values for Average Intensity (AvgInt) in the summer of 2012, 2014, 2016, 2017, 2018, and to a lesser degree 2013 and 2015. One may also notice that when Total Accumulation (TotAccum) was high and Duration was low, Maximum Intensity (MaxInt) was also high because it follows that higher intensity rainfall occurred when high rainfall totals were received in a short time-period. Time series of precipitation parameters. AvgInt and MaxInt refer to average and maximum precipitation intensity, respectively. TotAccum is the total depth of precipitation received during each weekly measurement period, and Duration is the total minutes during which precipitation was measured, for each measurement period. Columns delineate seasons (Su = summer, A = autumn, W = winter, and Sp = spring). Time series of precipitation parameters. AvgInt and MaxInt refer to average and maximum precipitation intensity, respectively. TotAccum is the total depth of precipitation received during each weekly measurement period, and Duration is the total minutes during which precipitation was measured, for each measurement period. Columns delineate seasons (Su = summer, A = autumn, W = winter, and Sp = spring).

Erosion
Mean erosion by measurement period (assessed using the average absolute change variables CAvg|Ch|, IAvg|Ch|, and SAvg|Ch|, where C, I, and S, refer to channels, interfluves, and sidewalls, Water 2020, 12, 925 6 of 15 respectively) was greater in winter and spring than the overall mean for all three geomorphic areas (Table 2). Notably, in winter months, CAvg|Ch| was 16.8 mm compared to 9.9 mm overall and SAvg|Ch| was 8.0 mm compared to 5.0 mm overall. Seasonal effects on interfluves were less pronounced, with IAvg|Ch| in winter at 4.8 mm compared to the overall mean of 3.5 mm. As with precipitation parameters, autumn was the season with the lowest mean erosion by measurement period for all geomorphic areas at 4.8 mm for channels, 3.5 mm for sidewalls, and 2.8 mm for interfluves. Seasonally, erosion variables show the most variability during winter months (Figure 3). Winter of 2016-2017 experienced less erosion than other years for all geomorphic areas, however, the study area received high rainfall accumulation during two weekly measurement periods.
Seasonally, erosion variables show the most variability during winter months (Figure 3). Winter of 2016-2017 experienced less erosion than other years for all geomorphic areas, however, the study area received high rainfall accumulation during two weekly measurement periods.

Statistical Modeling
Erosion variables were significantly correlated with total accumulation and duration parameters for all variables except interfluve erosion (IErosion) ( Table 3). Concordant with prior studies, erosion in channels was most strongly correlated with total accumulation (r = 0.467, r = 0.352, and r = −0.469 for CAvg|Ch|, CDep, and CErosion, respectively) and duration (r = 0.470, r = 0.367, and r = −0.447 for CAvg|Ch|, CDep, and CErosion, respectively). Note that all correlation coefficients for erosion variables (CErosion, IErosion, and SErosion) are negative because these variables are values below zero.

Statistical Modeling
Erosion variables were significantly correlated with total accumulation and duration parameters for all variables except interfluve erosion (IErosion) ( Table 3). Concordant with prior studies, erosion in channels was most strongly correlated with total accumulation (r = 0.467, r = 0.352, and r = −0.469 for CAvg|Ch|, CDep, and CErosion, respectively) and duration (r = 0.470, r = 0.367, and r = −0.447 for CAvg|Ch|, CDep, and CErosion, respectively). Note that all correlation coefficients for erosion variables (CErosion, IErosion, and SErosion) are negative because these variables are values below zero.
Spearman's correlation between the four precipitation parameters was compared to assess the potential for multicollinearity in statistical models, and total accumulation shows a very strong positive correlation with duration (r = 0.903) and a moderately strong positive correlation with average intensity (r = 0.591) and maximum intensity (r = 0.657). Likewise, average and maximum intensity were strongly and positively correlated (r = 0.794). Table 3. Spearman's correlation coefficients for erosion variables and precipitation parameters. C, channel; I, interfluve; S, sidewall. Only significant correlations are shown (* significant at α = 0.05, ** significant at α = 0.01).

Variable Name
Total Accumulation (mm) Before modeling erosion by season, OLS regression models were developed for the annual dataset (all measurement periods) using the four precipitation parameters from the current period, plus lagged variables for up to 11 prior periods (weeks). Table 4 summarizes output from models for each erosion variable in columns, with the variable name and R 2 value at the head of the column, and retained parameters marked by *. Retained parameters (independent variables) were those with statistically significant coefficients in each OLS model output. Nine models are represented in Table 4, one for each erosion variable. Model coefficients are not presented (only significance) here because the purpose of the modeling was to identify the precipitation parameters that were universally important, which was completed through frequency analyses. All model linear equations are, however, presented in Table A1 in Appendix A. Duration and total accumulation were the most important variables for channel erosion, while average intensity was important for erosion in interfluves and sidewalls. Also notable is the influence of antecedent precipitation at lags of up to 11 weeks for some variables. Seasonal OLS regression models clearly indicate the importance of precipitation intensity, which was, in prior studies, not retained in annual models of erosion ( Table 5). Note that seasonal models for IAvg|Ch| were omitted from Table 5 because only one viable model was generated, and its coefficient of determination was extremely low (R 2 = 0.064).
Interestingly, in summer and winter, average and maximum intensity were important explanatory parameters both during the current period, but also in prior periods. Precipitation intensity was not often retained in models of erosion during spring and autumn. It is also important to note that viable OLS regression models were generated for all erosion variables for summer and winter, with coefficients of determination ranging from R 2 = 0.245 to R 2 = 0.49 (except for IErosion in summer at R 2 = 0.131 and SDep in winter at R 2 = 0.087), suggesting that precipitation is an important driver for erosion in these months, no matter the metric used. Moreover, these results show that the character of the precipitation is an important driver for erosion; antecedent precipitation has an influence on erosion in the following weeks and months and it varies with season. Table 5. Parameters retained (indicated by *) in seasonal Ordinary Least Squares regression models of erosion variables (dependent variables) using lagged precipitation parameters (independent variables) Duration (min), Total Accumulation (TotAcc (mm)), and Average and Maximum Intensity (AvgInt and MaxInt, respectively (mm/min)). C, channel; I, interfluve; S, sidewall. Each column represents a separate model.

Erosion Variability
Variability exists in erosion statistics between the three geomorphic areas, such that channels had the highest variability and interfluves the lowest (Table 2), with sidewalls having intermediate variability. In particular, for both the overall annual dataset and for each seasonal partition, the mean and standard deviation were of similar magnitudes. Similar behavior was observed in a previous study in the same study area [13] and a study of gully erosion in the Karoo region of Africa [23]. Channels were dynamic and acted as both source and sink for sediment loads. Slugs of sediments gathered intermittently in the channel areas and were transported with channel flow following precipitation. Soil erosion was dominant in the gully sidewalls, however, the variability was moderate compared to channel erosion data, implying that sidewalls were less responsive with regard to erosion. In contrast, in the interfluve, the lesser amount of erosion and variability reflected the limited sediment yield, which may be due to the presence of vegetation that retarded erosion and lower gradient. Additionally, differences in soil cover thickness, soil types, moisture content, slope aspect and angle within the different geomorphic settings may explain the range of variability, however, that is beyond the scope of this paper and will be studied in the future.

Erosion-Precipitation Relationships
Seasonally, a comparison of erosion variables and precipitation parameters shows the same trend. Ordering seasonal precipitation parameters (Duration and TotAcc) and erosion variables from greatest to least, winter was greatest, followed by spring, summer and lastly, autumn. We see in Table 2 that winter months were the most dynamic, with the greatest mean erosion and the largest standard deviation of all seasons, and this pattern was consistent across channels, interfluves, and sidewalls for all erosion variables. This may be explained by the character of the winter precipitation: greater total accumulation and duration during these months associated with frontal precipitation events. Prior research has also demonstrated that freeze-thaw events are significant drivers of erosion in winter months at this site [14,19]. A similar pattern existed for spring, likely influenced by precipitation accumulation and duration as well as antecedent winter freeze-thaw activity [19]. Next, summer erosion and precipitation (Duration and TotAcc) ranked third, but interestingly, summer experienced the highest precipitation intensity of all seasons (both for AvgInt and MaxInt) ( Table 1). This reflected the dominance of convectional precipitation events in summer. Autumn experienced the minimum erosion and precipitation accumulation and duration, but greater maximum precipitation intensity than the annual average. This suggests that autumn precipitation events were short duration, high-intensity events that did not produce much precipitation depth and had little erosive power.
During winter 2016-2017, precipitation variables were near normal levels for the winter season, however, erosion for all geomorphic areas was very low (Figure 3). We examined temperature during this time period to determine whether the reduced freeze-thaw activity may have played a part, but, while winter 2016-2017 had less intense freeze-thaw activity than other winters during the study period, freeze-thaw events occurred. The timing of the greatest precipitation accumulation and duration was late autumn/early winter, and because these events were coincident, they indicate a period of low-intensity precipitation that may have encouraged more infiltration and less runoff, leading potentially to less erosion during this period. Lower hydrostatic pressure in unsaturated soils increases cohesion [24] which may be a significant factor associated with reduced erosion in late autumn and early winter of that year.
Average Intensity and Maximum Intensity of precipitation were very different, with approximately three orders of magnitude between the generally low average precipitation intensities and maximum intensity for the full dataset and each seasonal partition (Table 1). Future research at this site should assess the soil's infiltration capacity and explore different metrics that may better capture the relation between precipitation intensity and erosion. For example, measuring the rainfall duration when the rain rate exceeds the soil's infiltration capacity would generate a metric of the length of time during which there was a high probability of runoff generation.

Precipitation as a Driver for Erosion
Prior research at this site using 14 months of data found that Duration and TotAcc were the drivers for erosion, most strongly in channels. With six years of data, the present study confirmed the earlier result when erosion and precipitation data were lumped without regard for season. OLS regression models of annual erosion for the nine erosion variables, using the set of lagged precipitation parameters as independent variables, and overwhelmingly retained Duration parameters most frequently (24 times) ( Table 6). This means that overall nine OLS models of erosion outlined in Table 4, Duration and lagged Duration independent variables had significant coefficients 24 times. Despite the high correlation between TotAcc and erosion variables (Table 3), TotAcc was retained less frequently in the models (7 times) due to the high correlation between Duration and TotAcc (r = 0.903, p = 0.001) ( Table 3), indicating multicollinearity. Lagged intensity parameters were likewise retained fewer times; AvgInt parameters were retained 14 times, while MaxInt parameters were retained only 5 times. Therefore, using lumped annual data, Duration was the most important predictor of erosion, indicating that over the long term, prolonged precipitation is key. When erosion data were partitioned by geomorphic areas (Table 6), channel models overwhelmingly retained Duration most often. In contrast, sidewall and interfluve models retained Duration and AvgInt at approximately the same frequency (retained in 6 and 5 interfluve models and 8 sidewall models, respectively). This shows the importance of precipitation intensity as a driver for erosion in these two geomorphic areas. This may occur because interfluves and sidewalls may be more exposed to rain splash erosion, which is associated with higher intensity precipitation. Channels are not as steeply sloped as sidewalls and gully channel erosion is associated with the flow within the channel, which occurs after long-duration events that result in saturation-related runoff.
When erosion data were partitioned by season, the influence of precipitation intensity became apparent, especially during summer and to a lesser degree winter. This may be observed in Table 6, where MaxInt lagged parameters were retained 12 and 6 times in summer and winter erosion models, respectively, but only 0 and 3 times in spring and autumn models, respectively. This indicates that, while over the long term, Duration was the most important driver, during certain individual seasons intensity became important. This emphasizes the importance of the mechanics of convectional storms (summer) and frontal storms (winter) as an additional factor in seasonal erosion patterns. These patterns are also apparent when model results are partitioned by both season and geomorphic area ( Table 6).
Partitioning the data by season, therefore, produces additional knowledge that was not previously captured. We conclude that different drivers may be more effective agents of erosion in different seasons and, therefore, we recommend that studies of precipitation driven erosion should, wherever possible, partition data by season.

Conclusions
This study examined the effect of precipitation parameters on soil erosion through six years of high-resolution weekly monitoring in an Appalachian hillslope paying particular attention to seasonal effect. The long-term data provided an understanding of the seasonal pattern of soil erosion in a humid sub-tropical environment, which was not noticeable in other studies in the region using an annual dataset.
Different gully morphologies responded differently to long-term erosion. Channels were most active, showed a wide range of variability, and responded most dynamically, whereas the interfluves were least disturbed by erosion. Sidewalls were prone to erosion but were not as dynamic as channels.
To explore the reason behind varied gully erosion patterns in the different geomorphic settings, further studies are recommended to evaluate how erosion fluctuates with soil cover thickness, soil types, moisture contents, slope aspect, and slope angle.
Precipitation duration was the most important factor in initiating and continuing erosion year-round, yet seasonality played a significant role in the severity of gully erosion. Erosion was most pronounced in winter months, followed by spring, indicating the influence of high-intensity precipitation from frontal systems and repeated freeze-thaw cycles. Erosion in summer was driven by high-intensity precipitation from convectional storms. Soils in the study area were least prone to erosion during the moderate months of autumn. In channels, precipitation duration was the dominant driver for erosion due to runoff-related erosion, while in sidewalls and interfluves, intensity parameters were equally important as duration, likely related to rain splash erosion. This research shows that soil erosion is seasonally variable and an understanding of the seasonal pattern of soil erosion with respect to precipitation-related drivers improves the potential to achieve strategic conservation measures.  periods were retained most often, and these variables had the highest standardized coefficients compared to the intensity parameters (AvgInt and MaxInt) (standardized coefficients are not shown in the table). For interfluves, AvgInt and MaxInt were also retained in the models, and for the IDep and IErosion models, standardized coefficients for all retained variables were of similar magnitudes. For sidewalls, a similar pattern was generally noted, with retention of the intensity variables. For the SErosion model, Duration and TotAcc parameters had the largest standardized coefficients. Table A1. Regression equations for erosion variables (dependent variables) using lagged precipitation parameters (independent variables). Lagged variable names are appended with "LagN", where N indicates the number of measurement periods of antecedent lag. Duration_Lag1 indicates precipitation duration in prior measurement period (Lag of 1 period).