Quantifying the changes of soil surface microroughness due to rainfall-induced erosion on a smooth surface

This study examines the rainfall induced change in soil microroughness of a bare soil surface in agricultural landscapes. The focus is on the quantification of roughness length under the action of rainfall for initial microroughness length scales of 2 mm or less, defined here as initial smooth surface conditions. These conditions have not been extensively 15 examined in the literature as most studies have focused on initial disturbed surface conditions (bed surface conditions with initial length scales greater than 2 mm and varying between 5 – 50 mm). Three representative intensities namely 30 mm/h, 60 mm/h and 75 mm/h were applied over a smoothened bed surface at a field plot via a rainfall simulator. Soil surface microroughness measurements were obtained via a surface-profile laser scanner. Two indices were utilized to quantify soil surface microroughness, namely the Random Roughness (RR) index and the crossover length. Findings show a consistent 20 increase in roughness under the action of rainfall for initial microroughness length scales of 2 mm. This contradicts existing literature where a monotonic decay of roughness of soil surfaces with rainfall is recorded for disturbed surfaces. Analysis shows that on an average the RR and the crossover length post run increase by a multiple of 3.15 and 1.9, respectively from their corresponding values apriori the runs.


Introduction
Soil surface roughness influences many hydrologic processes such as flow partitioning between runoff and infiltration, flow unsteadiness, as well as soil mobilization and redeposition at scales ranging from a few millimeters to hillslope level (e.g.Huang and Bradford, 1990;Magunda et al., 1997;Zhang et al., 2014).
There are three distinct classes of microtopography surface roughness (Fig. 1a) reported in the literature for agricultural landscapes, each one of them depicting a representative length scale (Römkens and Wang, 1986;Potter, 1990).Following Oades and Waters (1991), the first class includes microrelief variations from individual soil grains to aggregates in the order Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2016Discuss., doi:10.5194/npg- -76, 2017 Manuscript under review for journal Nonlin.Processes Geophys.Published: 23 January 2017 c Author(s) 2017.CC-BY 3.0 License. of 0.053-2.0mm.The second class consists of variations due to soil clods ranging between 2-100 mm.The third class of soil surface roughness is systematic elevation differences due to tillage, referred to as oriented roughness (OR), ranging between 100-300 mm.
From the classes outlined above, the first two classes are the so called random roughness (RR), and constitute the main focus of the present research.Contrary to OR, which changes seasonally and during crop rotations, RR changes on an event base (Abaci and Papanicolaou, 2009).RR reflects the effects of rainfall action on the soil surface and inherently varies in space and time.As a result, RR affects key hydrologic processes at the soil scape and ultimately at the hillslope scale (e.g., overland flow, infiltration), by altering the depression storage and buffering or enhancing runoff and erosion processes (Gomez and Nearing, 2005;Chi et al., 2012).The most widely used statistical microrelief index for the evaluation of RR (Paz-Ferreiro et al., 2008) was initially proposed by Allmaras (1966) as the standard deviation of the log-transformed residual point elevation data.In this study, the RR index is calculated according to Currence and Lovely (1970) as the standard deviation of bed surface elevation data around the mean elevation, after slope and tillage effects have been removed in the individual height readings: where and ̅ are individual elevation height readings and their mean, respectively, and n is the total number of readings.
Characterization of RR remains a challenge.Most of the available studies are limited to soil surfaces where the length scale exceeds the upper microrelief length scale of 2 mm corresponding to the first class.These studies usually include bed surface conditions with initial length scales of 5 -50 mm (e.g., Zobeck and Onstad, 1987;Gilley and Finkner, 1991).In these studies, a monotonic decay of roughness due to precipitation action is predicted, since rainfall impact and runoff "smoothen" the rough edges of soil grains, aggregates and clods, especially in the absence of cover (Potter, 1990;Bertuzzi et al., 1990;Vázquez et al., 2008;Vermang et al., 2013).Few to none of existing studies, to the best of our knowledge, have examined however changes in RR for initial microroughness scales less than 2 mm.This condition remains prevalent usually between cyclic crop rotations in agricultural hillslopes, where the soil is "smoothened" from exposure to rainfall impact and runoff (Abaci and Papanicolaou, 2009).There is some experimental evidence that under such conditions, RR may actually increase under the action of rainfall (Huang and Bradford, 1992;Rosa et al., 2012;Zheng et al., 2014).We herein examine changes in RR under rainfall impact for initial microroughness less than 2 mm.In hillslope scale erosion models this condition is not represented adequately.In fact, most models assume a static roughness or at best a monotonic decay of RR, leading to errors on storage, runoff, and sediment routing (Liu and Singh, 2004).
This study expands our knowledge base for changes in RR, by investigating the interaction of rainfall with soil surfaces of microroughness length scales less than 2 mm.The main goal is therefore to examine the postulate observed in the literature (e.g., Zheng et al., 2014) that rainfall action on surfaces with initial microrelief less than 2 mm can lead to RR increase.An increase in RR under these scales would further imply that rainfall action cannot completely eliminate the roughness of a surface, so roughness residuals would always exist at the locations where raindrop detachment is prevalent.Hereafter, for shortness, tests with initial RR close to or less than 2 mm (first microroughness class) will be referred to as "smooth", whereas tests with initial RR greater than 2 mm will be referred to as "disturbed".
A key specific objective of this paper is to quantify changes in the microroughness of smooth bare soil surfaces under different rainfall intensities.We will first use the RR index because of its widespread use as a descriptor of soil microroughness both in soil microrelief studies and in the majority of existing landscape models.Additionally, we will supplement our RR analysis with the fractal properties of the soil surface by means of the crossover length derived from semivariogram analysis, taking advantage of its scale independence.The pros and cons of the two methods are discussed under bare soil with minimal microroughness.These findings are discussed with the aim to enhance current formulations by providing updates on microroughness during a rainfall event under bare soil conditions.

Experimental Conditions
This study was conducted on an experimental plot of the U.S. National Science Foundation Intensively Managed Landscapes Critical Zone Observatory in the headwaters of Clear Creek, IA (Fig. 1b).The soil series at the plot where the experiments were conducted is Tama (fine-silty, mixed, superactive, mesic Cumulic Endoaquoll) (http://criticalzone.org/iml/infrastructure/field-areas-iml/).The experimental plot was uniform in terms of downslope curvature, its gradient was 9% and the plot size was approximately 7 m long by 1.2 m wide.
The soil surface was prepared before each experiment by tamping using a plywood board to create a smoothened surface.
Rainfall was applied to the plots using Norton Ladder Multiple Intensity Rainfall Simulators designed by the USDA-ARS National Soil Erosion Research Laboratory, IN. Figure 2 shows the setup for all the experimental runs considered in the present study.For each test, three rainfall simulators were mounted in series over the experimental plot (Fig. 2a) and approximately 2.5 m atop the plot surface (Fig. 2b) in order to ensure that raindrop terminal velocity was reached.Water was continuously pumped from a water tank under controlled pressure, and uniform rainfall was applied through oscillating VeeJet nozzles which provided spherical drops with a median diameter of 2.25 mm and a terminal velocity of 6.8 m/s (see Table 1).The distribution of raindrop sizes generated by the rainfall simulators was calibrated using a disdrometer and followed a Marshall-Palmer distribution (Elhakeem and Papanicolaou, 2009), which is a widely accepted distribution for natural raindrop sizes in the U.S. Midwest where the study was performed.The calibration of the raindrop sizes was achieved by adjusting the pressure and swing frequency of the VeeJet nozzles.Special care was taken to minimize any potential biases with respect to natural rainfall, and, thus, render the rainfall simulation experiments scalable to other regions experiencing the same type of natural rainfall characteristics.
Surface elevations were obtained prior to and after the completion of the experiments via an instantaneous digital surfaceprofile laser scanner (Darboux and Huang, 2003), developed by the USDA-ARS National Soil Erosion Research Laboratory, IN (Fig. 3a).Laser scanner measurements a priori the runs ensured that the overall microrelief was less than 2 mm.After the assembly and their movement on the carriage is controlled from the desktop using a computer program that regulates the travel distance based on a user-specified distance (Fig. 3a).Information captured by the camera is recorded with an attached computer.The information from each scan is converted into a set of (x,y,z) coordinates using a calibration file and specific software for data transformation.The set of (x,y,z) coordinates obtained for each experiment are imported into ArcGIS 10.3.1 in order to create the corresponding Digital Elevation Models (DEMs) through inverse distance weighting interpolation and thereby visualize or analyze the surfaces (Fig. 3b).The resulting DEMs have a horizontal resolution of 1 mm and an accuracy of 0.5 mm in the vertical.
Three experimental tests of varying rainfall intensity were conducted on the experimental plot.Experiment 1 had a rainfall intensity of 30 mm/h, Experiment 2 an intensity of 60 mm/h, and Experiment 3 an intensity of 75 mm/h (Table 1).These simulated intensities represent typical storms observed in the region of South Amana where the plot is located.The initial soil water content was measured with Decagon soil moisture sensors and was the same for each experiment and approximately equal to 35% at the whole plot.Each experiment was run for nearly 5 hours, sufficiently long to reach steady state runoff and sediment, as confirmed by weir readings and discrete samples taken at the outlet of the plot.Notice that the initial microroughness length scale in Experiment 1 (1.17 mm) was greater than that of Experiment 2 (0.42 mm) and Experiment 3 (0.32 mm)see Table 2.This is attributed to the different timing of the experiment runs with respect to tillage.Experiment 1 was performed in early August, soon after harvest, so the soil surface had recently been disturbed.
However, for Experiments 2 and 3 which were performed in late September, the soil presented less surface disturbance due to the cumulative action of rainfall within that period (Papanicolaou et al., 2015).Therefore, despite tamping with plywood remnants of tillage effects remained in Experiment 1 yielding different initial microroughness length scales than Experiments 2 and 3.All cases, however, exhibited initial microroughness length less than 2 mm corresponding to smooth surface bed conditions as confirmed with the laser scanner.Dry soil bulk density was 1.25 kg/m 3 for Experiment 1, and about 6% higher for Experiments 2 and 3 due to self-weigh consolidation of soil.
The left part of Fig. 4 provides an example of the experimental plot at pre-rainfall and post-rainfall conditions.Since the focus of this research is only on plot regions where raindrop detachment is dominant over runoff, we are using the scanned profiles that correspond only to these upslope locations, which are shown at the right part of Fig. 4. By definition, no rill formation ever took place in these regions.For scanned profiles within the Region of Interest (ROI) (i.e., a selected 200 mm x 200 mm window size), we extract the data for further statistical and geostatistical analyses by utilizing the public domain R package (https://www.r-project.org/).The geostatistics (gstat) and spatial analysis (sp) libraries are imported to create sample semivariograms.

Soil Surface Roughness Quantification
Due to its commonality found in the literature and its widespread use in various landscape models as a descriptor of microroughness, the RR index calculated from Eq. ( 1) is used in this study as the first method to describe soil surface roughness.The RR index, however, requires that there is no spatial correlation between the surface elevations (Huang and Bradford, 1992).If correlation exists within a certain spatial scale, RR will likely change with the changing window size of observed data (Paz-Ferreiro et al., 2008) and may be dependent on the resolution of the measurement device (Huang and Bradford, 1992).Thus, alternative scale-independent indices from methods that consider spatial correlation have been considered in order to address this issue.These methods include semivariogram analysis (Vázquez et al., 2005;Oleschko et al., 2008;Rosa et al. 2012;Vermang et al., 2013), fractal models based on Fractional Brownian Motion (Burrough, 1983a;Vázquez et al., 2005;Papanicolaou et al., 2012;Vermang et al., 2013), multifractal analysis (Lovejoy and Schertzer, 2007;Vázquez et al., 2008), Markov-Gaussian Model (Vermang et al., 2013), and two-dimensional Fourier Transform (Cheng et al., 2012), among others.
The crossover length derived from semivariogram analysis is an index that is commonly used in most recent soil microrelief studies to describe surface microroughness, with the advantage of its quantification being scale independent through the consideration of the spatial correlation between surface elevations (Vázquez et al., 2007;Paz-Ferreiro et al., 2008;Tarquis et al., 2011).The semivariogram is a useful geostatistical tool developed to depict the spatial autocorrelation of data.It is calculated from the following equation: where ( ) is the semivariance, is the lag-distance between data points, ( ) is the elevation height value at location and ( ) is the total number of pairs separated by lag-distance considered in the calculation.The semivariogram is the plot of the semivariance with respect to the lag-distance.
Key fractal indices for describing soil surface roughness can be derived from the semivariogram.Assuming a fractional Brownian motion model for describing soil surface roughness, as proposed in the pioneering work of Mandelbrot and van Ness (1968), the following expression for ( ) that incorporates the generalized Hurst exponent, is obtained (Huang and Bradford, 1992;Vázquez et al., 2007;Paz-Ferreiro et al., 2008): where H is a measure of the degree of correlation between the surface elevations at lag distance h and l is the crossover length and .The crossover length is a measure of the vertical variability of soil surface roughness at the particular scale where the fractal dimension is estimated, hence greater roughness is associated with larger crossover length values and vice versa (Huang and Bradford, 1992).
Given the semivariogram plot from Eq. ( 2), H and l can be extracted by fitting a power law relationship in the form of to the semivariance-lag distance data.The B regression variable gives the generalized Hurst exponent value, whereas the A regression variable yields the crossover length.In this study, a 3 rd order polynomial surface is optimally fitted within the predefined ROI, and RR is calculated as the standard deviation of the residuals from the aforementioned best fit surface based on Eq. ( 1).The semivariograms for the experiments confirm the suitability of the 200 mm square window in capturing surface elevations at a scale where no spatial correlation exists, as explained in the results below.

Results
This section presents the results from the rainfall experiments in the form of ratios of the roughness indices and it is organized as follows.First, the RR ratio, defined as the ratio of the RR index post-rainfall over the RR index prior to the rainfall (RR post /RR pre ), is calculated for each experiment and compared to existing literature.Next, semivariograms are plotted under pre-and post-rainfall conditions at the ROI to assess the spatial correlation of surface elevations and the associated crossover lengths are calculated.Finally, the crossover length ratio, defined as the ratio of the crossover length post-rainfall over the crossover length prior to the rainfall (l post /l pre ), is calculated for each experiment, compared to existing literature, and assessed as a descriptor of change in microroughness along with the RR ratio.The results are presented in the form of ratios to negate the effects of the differences found in initial microrelief amongst the three runs, thereby comparing changes in relative terms.

Changes in the RR index
Based on visual inspection of the DEMs on the right hand side of Fig. 4, it is evident that roughness at the upslope regions increases with rainfall.Figure 5 provides quantitative information and shows the RR ratio, i.e., RR post /RR pre , which describes the change in roughness estimated within the ROIs with respect to the initial microroughness conditions.From Fig. 5 it can be seen that the resulting RR index is 1.34, 3.55 and 4.56 times higher than its initial value when rainfall is applied for Experiments 1, 2 and 3, respectively (also see Table 2).
Table 2 also summarizes the effects of rainfall on roughness by means of the RR index for the present study along with other representative studies in existing literature, capturing the range of smooth and disturbed initial microroughness conditions found in the literature.The last column shows the RR ratio.From Table 2, two main observations can be made.First, for the smooth surface initial condition studies a consistent increase in roughness with precipitation occurs.Second, the present study indicates that the RR index becomes higher with higher rainfall intensity (Fig. 5, Table 2).
The results outlined above for the use of the RR index as a descriptor of change in microroughness have been based on the assumption that there is no statistically significant spatial correlation in elevation readings between neighboring locations at the ROI, so they are valid only under this assumption.The following subsection outlines and discusses the results of the semivariogram analysis and the crossover length in order to confirm this assumption and compare with the RR index method.

Changes in the crossover length
Semivariograms were obtained from geostatistical analysis and plotted at four different angles -0°, 45°, 90°, and 135°-with respect to the downslope direction in order to test isotropy.The root mean square error of the semivariogram values at prerainfall and post-rainfall conditions was found to be less than 10%.Therefore, there were no statistically significant differences in the semivariograms, which implies that the roughness patterns for smooth initial conditions at the ROI remain isotropic after a rainfall event and can thus be examined with one representative semivariogram.
The representative sample semivariogram in the downslope direction calculated at the ROI for each experiment are presented in Fig. 6.The vertical dashed lines designate the lag distances above which the spatial autocorrelation of the elevations is not statistically significant.These lag distances are approximately 10 mm, so the selected 200 mm window size of the ROI is almost 20 times greater than the spatial autocorrelation range, which is considered sufficient to assume no spatial autocorrelation at the scale examined in this study.Hence, RR is independent of the window size of observed data considered.
Moreover, the semivariogram sill (which is defined as the near-constant value of semivariance at large lag distances where the semivariogram levels outsee horizontal dashed lines in Fig. 6) is directly related to the RR of the soil surface since it quantifies the mean variance in elevations at distances where spatial correlation is statistically negligible (Vázquez et al., 2005).It is seen from Fig. 6 that the post-rainfall sills are greater than their corresponding pre-rainfall value for all three intensities.Also, the difference in sills between pre-and post-rainfall conditions for the 30 mm/h precipitation intensity is much lower than those of the 60 mm/h and 75 mm/h precipitation intensities.These observations are in accordance with the results noted earlier for the RR ratio (see Fig. 5, Table 2).Therefore, similar conclusions drawn from the RR ratio graphs can be drawn from the semivariograms.
Table 3 summarizes the impact of rainfall on the crossover length l at the ROI for our present study along with other representative studies in existing literature.As seen from the table, the final roughness state of the bed surface after the application of rainfall (for all the studies considered herein), tends to have a crossover length in the range of 0.2 -4 mm.
However, the existing studies with initial disturbed surface conditions report a decrease in the crossover length after rainfall, contrary to the observed increase in the crossover length reported in our study for initial smooth conditions (e.g., Vázquez et al., 2007;Paz-Ferreiro et al., 2008;Vermang et al., 2013).Fig. 7 shows the crossover length ratio at the ROI with respect to rainfall intensity for our study.Crossover length ratios greater than unity reflect an increase of soil surface roughness with rainfall.Similar to the RR ratio, the crossover length ratio is greater at the high precipitation intensity cases (Experiments 2 and 3) than at the low precipitation intensity case (Experiment 1)also reflected in Table 3. Overall, these results suggest that the crossover length and RR index, for scales where no spatial correlation exists, can be used interchangeably in order to characterize rainfall induced changes in microroughness.

Conclusions and Discussion
We present unique, novel rainfall simulation experiments in this study that mimic natural rainfall conditions at the plot scale and capture microroughness evolution under bare surface conditions with microroughness in the order of 2 mm.Analysis of soil surface roughness in the region where raindrop detachment dominates and under initial smooth surface preconditions for three rainfall intensities shows a consistent increase in the RR index and crossover length, which are confirmed as reliable descriptors of microroughness.This increase contrasts the findings of most available literature and is attributed to the initial roughness magnitude examined in this study; that is, smooth surface preconditionsdefined herein as conditions of microroughness length scale close to 2 mm or less (first microroughness class)as opposed to the existing literature that mostly deals with initial minimal roughness magnitudes in the range of 5-50 mm.Our findings indicate the existence of a characteristic roughness scale in the magnitude of 2 mm below which RR is expected to increase and above which RR is expected to decrease due to the action of rainfall.Another significant outcome of this study is the fact that the mere action of rainfall cannot completely smoothen out a bed soil surface, thereby a roughness residual will always remain at the locations where the action of runoff is low or absent.Roughness residuals further infer depression storage residuals, which can significantly alter the flow partitioning between infiltration and runoff at the hillslope scale especially at the onset of an event (Onstad, 1984).This finding benefits a number of disciplines, from environmental modeling and research, to agricultural management, by providing a better understanding of the highly dynamic phenomenon of soil surface microroughness evolution.For instance, current modeling tools of soil surface processes are likely to predict a total decay of soil surface roughness after subsequent rainfall events (e.g., Potter, 1990) that will be translated to the absence of depression storage and are prone to overpredicting runoff (Liu and Singh, 2004).This may have implications on belowground processes within the soil column such as porous network formation and water nutrient diffusion (Hunt and Ghanbarian, 2016).Finally, our study demonstrates that soil surface roughness, which may be quantified by means of the RR index or crossover length, evolves in response to and in concert with rainfall, indicating a feedback where hydrologic response and surface evolution are interrelated.Our study's results further suggest that the response of soil surface roughness is nonlinear with respect to the hydrodynamic forcing and the related processes, depending on both the initial roughness condition and the applied rainfall energy and affecting ponding, flow pathways, as well as the associated surface heterogeneity.
On an annual basis, Abaci and Papanicolaou (2009) and Abban et al. (2016) highlight the importance of the seasonal variation of land cover on sediment output in agricultural Intensively Managed Landscapes (IMLs), indicating that during certain periods, the combination of high magnitude events and bare soil will severely increase erosion.This point is of high relevance here given the soil surface in agricultural IMLs is bare 75% of the time during the calendar year.Models simulating these periods are likely to be sensitive to the treatment (and definition) of the soil surface microroughness, and thus, require an adequate determination of the soil surface roughness length scales for accurately modeling the hydrologic response of hillslopes.They must be able to adequately capture the increasing and decreasing trends in soil surface microroughness to accurately predict the landscape response to precipitation.From this standpoint, our study provides a Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2016Discuss., doi:10.5194/npg--76, 2017     Manuscript under review for journal Nonlin.Processes Geophys.Published: 23 January 2017 c Author(s) 2017.CC-BY 3.0 License.
Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2016Discuss., doi:10.5194/npg--76, 2017     Manuscript under review for journal Nonlin.Processes Geophys.Published: 23 January 2017 c Author(s) 2017.CC-BY 3.0 License.runs the measurements provided changes in microroughness under the action of rainfall.Horizontal and vertical accuracies of the laser are 0.5 mm.Points were measured every 1 mm.The system consists of two laser diodes mounted 40 cm apart to project a laser plane over the targeted surface.The beam is captured by an 8-bit, high-resolution progressive scan CCD camera with 1030 rows x 1300 columns and a 9 mm lens.The camera and lasers are mounted on a 5 m long carriage Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2016Discuss., doi:10.5194/npg--76, 2017     Manuscript under review for journal Nonlin.Processes Geophys.Published: 23 January 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 1 :
Figure 1: (a) Definition sketch of soil surface roughness types.(b) Experimental plot where the experiments were conducted.The rainfall simulator is placed above the bare soil surface and a base made of wood is put into place to facilitate the movement of the surface-profile laser scanner.

Figure 2 :
Figure 2: Sketch of the setup for all of the experimental tests considered in this study.(a) The rainfall simulators are mounted in series and a pump is utilized to provide them with water from a water tank throughout the experiments.(b) The rainfall simulators are placed and adjusted at a height of 2.5 m above the experimental plot surface to ensure drop terminal velocity is reached.5

Figure 3 :
Figure 3: (a) Instantaneous digital surface-profile laser scanner used in the experimental runs and laser beam projected on the soil surface.(b) Cloud of (x,y,z) data acquired from the laser scanner for an experimental test along with the associated 3D representation of the soil surface microrelief through interpolation.5

Figure 4 :
Figure 4: Left: Part of the experimental plot of this study under pre-and post-rainfall conditions for an experimental test.The dashed boxes indicate the extent of the ROI.Right: Scanned profiles extracted from the laser-scanned areas of the three experimental tests considered, under both pre-and post-rainfall 5

Figure 5 :
Figure 5: Random Roughness (RR) Ratio at the ROI for the three experimental tests considered herein.

Figure 6 :
Figure 6: Spatial semivariograms at the ROI for the three experimental tests considered herein, under pre-5

Figure 7 :
Figure 7: Crossover length (l) ratio at the ROI for the three experimental tests considered herein.

Table 2 : Summary of the rainfall induced change of the RR index in the experimental tests of this study, as well as in experiments reported in the literature. Smooth conditions refer to initial microroughness close to or less than 2 mm and disturbed conditions refer to initial microroughness greater than 2 mm.
Nonlin.Processes Geophys.Discuss., doi:10.5194/npg-2016Discuss., doi:10.5194/npg--76,2017Manuscript under review for journal Nonlin.Processes Geophys.Published: 23 January 2017 c Author(s) 2017.CC-BY 3.0 License.