Nonlinear Processes in Geophysics Quantification of tillage , plant cover , and cumulative rainfall effects on soil surface microrelief by statistical , geostatistical and fractal indices

Changes in soil surface microrelief with cumulative rainfall under different tillage systems and crop cover conditions were investigated in southern Brazil. Surface cover was none (fallow) or the crop succession maize followed by oats. Tillage treatments were: 1) conventional tillage on bare soil ( BS), 2) conventional tillage ( CT ), 3) minimum tillage (MT ) and 4) no tillage ( NT ) under maize and oats. Measurements were taken with a manual relief meter on small rectangular grids of 0.234 and 0.156 m 2, throughout growing season of maize and oats, respectively. Each data set consisted of 200 point height readings, the size of the smallest cells being 3 ×5 cm during maize and 2 ×5 cm during oats growth periods. Random Roughness ( RR), Limiting Difference (LD), Limiting Slope (LS) and two fractal parameters, fractal dimension ( D) and crossover length (l) were estimated from the measured microtopographic data sets. Indices describing the vertical component of soil roughness such as RR, LD andl generally decreased with cumulative rain in theBS treatment, left fallow, and in the CT and MT treatments under maize and oats canopy. However, these indices were not substantially affected by cumulative rain in theNT treatment, whose surface was protected with previous crop residues. Roughness decay from initial values was larger in theBS treatment than inCT andMT treatments. Moreover, roughness decay generally tended to be faster under maize than under oats. The RR and LD indices decreased quadratically, while the l index decreased exponentially in the tilled, BS, CT andMT treatments. Crossover length was sensitive to differences in soil roughness condiCorrespondence to: E. Vidal Vázquez (evidal@udc.es) tions allowing a description of microrelief decay due to rainfall in the tilled treatments, although better correlations between cumulative rainfall and the most commonly used indicesRR andLD were obtained. At the studied scale, parametersl andD have been found to be useful in interpreting the configuration properties of the soil surface microrelief.


Introduction
The Earth's topography encompasses a considerable range of scales.The horizontal scale of lithosphere relief varies from sizes smaller than millimeters to sizes as large as the planet with a perimeter in the order of 40 000 km (Huang, 1998;Lovejoy and Schertzer, 2007).The vertical scale of lithosphere relief reaches some 18 km considering fluctuations between oceanic bed and continental mountains.Therefore, topography variability on our planet is higher than factors of over 10 10 and 10 7 in the horizontal and vertical scales, respectively.At the planetary scale, topography remains relatively unchanged except in active volcanic and tectonic areas (Lovejoy and Schertzer, 2007).
Soil microrelief is referred to as the small scale topographic variation across a cultivated field (Allmaras et al., 1966;Römkens and Wang, 1986;Huang and Bradford, 1992).Conversely, soil roughness has been defined as "a measure of variation in surface microtopography" (Huang, 1998).For assessing soil surface microtopography features, point elevation readings are currently taken with millimeter to centimeter resolution within a meter scale area.This also provides an operative definition of soil surface microtopography (Huang, 1998;Vidal Vázquez et al., 2005).
Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.576 J. Paz-Ferreiro et al.: Microrelief fractal parameters Both landscape topography and soil surface microrelief are the result of many competing influences.However, unlike topography at the landscape scale, the microrelief can change rapidly in agricultural fields in response to tillage operations, erosion, deposition and other less dominant factors.Moreover, the roughness of the soil surface results from different elements such as individual grains, aggregates, clods, tillage marks and landscape features that each contribute at various scale lengths.Conversely, soil roughness can be split into several components.First, surface variations due to the spatial distribution of individual grains, aggregates and clods, define random roughness.Second, systematic differences in elevation caused by tillage implements such as furrows are referred to as oriented roughness; if concentrated erosion occurs, rills and gullies also contribute to oriented roughness.Third, at the field and landscape scales higher order roughness created by the general topographic slope with its concavities and convexities may be recognized; this type of roughness is the most interesting to sciences concerned with topographic or macrorelief scales (Römkens and Wang, 1986).
Soil surface roughness affects many transfer processes on and across the soil-atmosphere boundary, for example, infiltration, runoff, soil detachment by water and wind, gas exchange and evaporation.The temporal water storage in microrelief depressions controls overland flow generation and runoff pathways (Hairsine et al., 1992;Hansen et al., 1999;Favis-Mortlock et al., 2000).The quantification of roughness related parameters such as depression storage is considered of great relevance and practical importance, because it can be used to enhance both water conservation and soil conservation.
Soil surface roughness can also act as an erosivity factor at smaller scales within a microplot, affecting the soil response to raindrop impact and other erosive forces, and as an erodibility factor at larger scales of the sizes of a plot, a field or a hillslope, affecting the origin and spatial organization of the drainage pattern (Merril et al., 2001).This demonstrates the interaction between the length scale of a physical process on the soil surface and the scale of roughness components.
Because the random roughness condition is initiated by the variable disposition of aggregates and clods on the soil surface roughness, it is associated to disordered microrelief (Huang, 1998;Vidal Vázquez et al., 2005).Therefore, the term "random" only describes non-orderly distribution of structural unities on the soil surfaces and should not be misleading, since random roughness is spatially correlated at short distances (Linden and van Doren, 1986;Huang and Bradford, 1992;Helming et al., 1998;Eltz and Norton, 1997;Miranda, 2000;Zribi et al., 2000;Vidal Vázquez et al., 2005, 2006, 2007).This correlation seems to be a usual property of soil surface random roughness irrespective of soil type and tillage condition.The range of spatial dependence is of the same order of magnitude as the largest structural unities, aggregates and clods, on the soil surface.Above distances larger than the main structural unities, correlation tends to disappear.
Oriented roughness can be easily quantified by deterministic models, because it depends on slope and/or on a specific tillage tool (Huang, 1998).In contrast, random roughness analysis has been more challenging, as it should include a characterization of the spatial structure of microrelief elements with various sizes.Soil surface roughness, however, has been most frequently described by the random roughness or the root mean square roughness index (Allmaras et al., 1966).This index is equivalent to the standard error of point elevations relative to a reference plane and represents a statistical measure of vertical topographic fluctuation, implicitly assuming that there is no spatial dependence in surface roughness.But its main shortcoming is that two soil surfaces with identical RR values may have different topographies (Römkens and Wang, 1986;Merril et al., 2001).
Furthermore, soil surface roughness, measured as RR, increases as the area extent of point elevation measurements is increased.Consequently, roughness should be expressed as a scale-dependent function, not merely as a statistical index (Huang, 1998).Different approaches have been proposed to account for the spatial scale.Based on the first order semivariogram, i.e., elevation difference measure, Linden and van Doren (1986) defined two indices: limiting difference (LD) and limiting slope (LS), describing the vertical component and the separation length scale of soil surface roughness, respectively.Correlation and spectral density analysis have been used to quantify roughness and to examine roughness changes as a function of rain and tillage operations (Currence and Lovely, 1970;Dexter, 1977).
Fractal models have been used to describe not only the Earth's topography (e.g.Mandelbrot, 1983;Lovejoy and Shertzer, 2007) but also soil microrelief (Armstrong, 1986;Huang and Bradford, 1992;Eltz and Norton, 1997;Huang, 1998;Miranda, 2000;Vidal Vázquez et al., 2005, 2006, 2007).The term "fractal" is often taken to be synonymous with "scale-invariance", the well known property that many geologic and soil characteristics look the same at all scales.This is referred to as self-similarity.In a self-similar fractal the factor that characterizes the invariance is independent of the direction.However, a fractal is identified as selfaffine when different scaling ratios are found for each independent direction.Early studies of Mandelbrot (1985) and Voss (1985) showed that the height of topography along a linear track can only be self-affine and not self-similar.By extension, topographic surfaces are also examples of self-affine fractals (Vidal Vázquez et al., 2005;Lovejoy and Shertzer, 2007).
Although few studies have addressed the quantification of soil surface roughness by fractal models, different techniques have been proposed for estimating fractal dimension, D, from point elevation data sets.Principally, two different categories of fractal models have been applied for assessing soil microrelief: variational and non-variational models (Vidal Vázquez et al., 2005).It should be emphasized that non-variational techniques such as tortuosity evaluation (Bertuzzi et al., 1990) and Richardson number (Gallart and Pardini, 1996;Pardini and Gallart, 1998;Pardini, 2003) are only appropriate for self-similar surfaces and therefore they should not be applied for estimating fractal parameters from point elevation data sets, which are self-affine.In contrast, variational models allow a description of soil surface microrelief from self-affine data sets.Because variational and non-variational models are not equivalent, results from both of these types of methods frequently do not agree (Miranda, 2000;Vidal Vázquez et al., 2005) and cannot be directly compared.
Two parameters are required to describe the scale dependent roughness function by a variational fractal model, namely fractal dimension, D, and crossover length, l.Fractal parameters, D and l, account for the multiscale effects and for the fluctuations of local vertical statistics, respectively (Huang and Bradford, 1992;Eltz and Norton, 1997;Huang, 1998;Miranda, 2000;Vidal Vázquez et al., 2005, 2007).
Soil surface roughness controls microrelief depression pattern and is clearly relevant to estimate depression storage capacity (Huang and Bradford, 1990;Kamphorst et al., 2000;Darboux et al., 2005).Therefore, the scale dependence of soil microrelief has implications in modeling this parameter.However, most runoff and erosion processes are not affected by the short-distance correlation so they could be quantified by a single parameter, without taking into account these scale dependent function.In other words, fractal-based roughness parameters are not necessarily more important than other indicators, based on statistics or on geostatistics.
This paper aims to clarify the relevance of the fractal approach for soil surface roughness evaluation.It extends previous studies on data sets acquired with low technology devices, characterizing large scale surface features.The main focus is on soil management effects over initial roughness created during tillage operations and on roughness destruction and decay by cumulative rainfall under a crop succession.

Site, soil and tillage operations
This study was conducted at the agricultural research station of the University of the State of Santa Catarina (UDESC) in Lages, Santa Catarina State, Brazil, latitude 27 • 49 S, longitude 50 • 20 W and mean altitude 937 m a.s.l.The climate was classified as Cfb type according to Köppen, with an average yearly rainfall and erosivity of about 1600 mm and 6000 MJ mm ha −1 h −1 , respectively (Bertol et al., 2002a).The experimental area has been used for water erosion studies under natural rainfall since November 1988.The data sets for our study were acquired between May 2003 and October 2004.
The studied soil with a moderate A horizon and a silty-clay sedimentary bedrock, was classified as an Inceptisol (Soil Survey Staff, 1993).The topsoil (0-20 cm depth) is clay textured, with 142 g kg −1 sand, 437 g kg −1 silt and 421 g kg −1 clay contents.Bulk density in this layer was 1.28 kg dm −3 .Soil erodibility was 0.0115 t ha h ha −1 MJ −1 mm −1 (Bertol et al., 2002b) and the mean slope of the experimental plots was 0.10 m/m.
Field experiments were conducted during the growth period of two successive crops, maize and oats, using 3.5 m×22.1 m erosion plots.Tillage treatments were as follows: 1) ploughing and two times harrowing on bare soil, as a control treatment (BS), 2) conventional tillage which consisted also of ploughing followed by two successive harrowing (CT ), 3) minimum tillage which consisted of chisel ploughing and harrowing (MT ) and 4) no-tillage over previous crop residues (N T ).Tillage systems were in a randomized design, without replications.Treatments CT , MT and N T were cropped first to maize (Zea mays) and after to oats (Avena strigosa), both crops being a part of a multiple rotation system.Cumulative rain amounted to 229 mm and 350 mm in the growth periods of maize and oats, respectively.A detailed description of tillage operations and experimental setup was presented elsewhere (Bertol et al., 2002(Bertol et al., , 2006)).
Operations with different tillage implements produced soil surfaces visibly different between treatments.In the conventionally tilled treatments (BS and CT ) as well as in the minimum tillage treatment (MT ) a range of aggregates and clod sizes was observed.These structural elements were more or less evenly distributed in the BS and CT treatments.However, in MT treatment rougher areas of disturbed soil due to chisel ploughing were distinct from undisturbed smoother areas between the chisel rows.Furthermore, CT and MT treatments were partly covered by previous crop residues that contributed to microrelief, because they were not totally incorporated into the soil by ploughing.Therefore, partial residue soil cover was larger in MT than in CT treatment.Under N T soil surface microrelief consisted of undulated land without tillage marks, fully protected by free crop residues, covering the entire soil surface.

Field data sets and trend removal
Soil surface microrelief was measured five times in the maize growth season and four times in the oats growth season.The first measurement was made just at sowing time; then successive readings were taken at about 15 days apart throughout each crop growth period.
Soil microrelief was measured with a pin meter, which held 40 calibrated needles.The interval between surface elevation readings made in the maize growth period was set at 3 cm, thus extending over a length of 117 cm.However, 578 J. Paz-Ferreiro et al.: Microrelief fractal parameters measurements in the oats season were performed with needles spaced 2 cm apart, along 78 cm transects.In both cases 5 profiles were taken per plot, 5 cm apart from each other.This gave 200 readings at any location.Therefore, the sampling scheme was a 3 cm×5 cm rectangular grid when the soil was cropped to maize and it was a 2 cm×5 cm grid in the oats season.The area of the roughness experimental plot was 0.234 and 0.156 m 2 in the maize and oats crops, respectively.For each transect, the pin position was registered photographically and later digitized to record readings of 40 calibrated needles, as described by Wagner and Yiming Yu (1991).
Manual pinmeters are destructive devices.Thus, different microplots were used for surface roughness measurements at increasing amounts of cumulative rainfall in successive dates.Within each tillage treatment, experimental plots were located as close as possible to minimize the effect of spatial variability between them.
Microrelief data sets were corrected for slope and tillage marks.Correction for slope was obtained using the plane of best fit to the 200 point elevation readings of each plot.Additionally, non-directional random roughness surfaces were obtained by removing row and column effects, as proposed by Currence and Lovely (1970).Therefore, all the indices in this study were estimated for data sets representing the random roughness condition.
The configuration of soil topography was simple described by a set of points of known x-, y-and z-coordinates.The elevation values given as a function of the horizontal coordinate system provide a numerical representation of the surface and constitute a digital elevation model (DEM).From each experimental data set of soil surface microtopography a DEM was obtained after trend removal, representing the random roughness condition.

Statistical and geostatistical indices
Besides quantifying of fractal parameters, in this study three traditional roughness indicators were assessed: a statistical index, random roughness (RR), and two geostatistical indices estimated from the first order semivariogram of point elevation differences, referred to as "limiting difference" (LD) and "limiting slope" (LS).
In accordance with Currence and Lovely (1970) and Kamphorst et al. (2000), RR was estimated simply as the standard deviation of height readings after correction for slope and tillage marks.Therefore, in this study RR was assessed without a log-transformation of the residual point elevation data and without removing extreme values, as initially proposed by Allmaras (1966).Random Roughness was calculated as: where Zi and Z are point elevation and mean point elevation, respectively, and n is the number of experimental points.
LD and LS quantification is based on the first-order variogram of mean absolute differences in elevation ( Z h ) separated by a distance vector, h.Following Linden and van Doren (1986), the first step in estimating LD and LS consists in the computation of the mean absolute elevation difference, as: where Z i and Z i+h are elevations at positions located a horizontal distance h apart and n is the number of data points.Therefore, for a two-dimensional data set, Z i+h corresponds to point heights located in a disk of radii h around point i.
Then, parameters, a and b, of the linear relationship between 1/ Z h and 1/ h are estimated: Parameters LD and LS are quantified, respectively, as: The limiting difference, LD, is mathematically interpreted as the value of Z when h approaches large spatial intervals.The limiting slope, LS, is the slope or Z/ h, when h approaches zero.This analysis is essentially independent of the minimum spacing interval because of the linear nature of Eq. (3).
In estimating roughness indices based on geostatistical concepts (LD and LS), surface microrelief was assumed to be statistically homogeneous, which implies that statistical properties do not depend on the position, but only depend on the spatial separation, h.

Fractal dimension and crossover length
A fractal analysis was performed on the soil microrelief data sets using a variational method, thus assuming a self-affine model for microrelief description.This method is based on a fractional Brownian motion (fBm) model to calculate the Hurst exponent, H , from which the fractal parameters D and l are obtained.
The semivariance function or semivariogram, γ (h), was selected as structural o scale-dependent function, because of the rectangular shape of the sampling grid (Miranda, 2000).Under these conditions, the semivariogram is more accurate and gives more consistent results than other structural functions such as the root-mean square (RMS) algorithm (Malinverno, 1990;Miranda, 2000).The semivariance can be estimated from sampling data as: where γ ×(h) is known as the experimental semivariogram, n is the number of pairs of sample points of observations of the values of the studied attribute, Z, between any two places x and x+h, separated by a distance h, Z(x) is the elevation at location x.Data sets were checked for the intrinsic hypothesis of the regional variable theory.After trend removal the two conditions defining the requirements for the intrinsic hypothesis, i.e., stationarity of differences and variance of differences were verified.
The use of the γ (h) function as the roughness measure gives a direct relationship between an elevation difference term and the separation length scale, h.This is an important step in fractal analysis as the estimated values of fractal parameters, D and l from soil topography surfaces or transects, may show some bias depending not only on the assumptions made in formulating the fractal model but also on the algorithm which is used (Miranda, 2000;Vidal Vázquez et al., 2005, 2006).Burrough (1983 a,b) first used the γ (h) structural function in soil science.The semivariance function for estimating fractal dimension of soil height tracks was introduced by Armstrong (1986) and later on applied to various surface types (Carr and Benzer, 1991;Davis and Hall, 1999), including agricultural soils (Huang and Bradford, 1992;Miranda, 2000;Vidal Vázquez et al., 2005, 2006, 2007).
The fractal model used for soil surface roughness quantification was the fractal Brownian motion model (fBm).The elevation difference of a fBm is given by: where the exponent of the incremental function, H , is the Hurst exponent.The power model, which describes a selfsimilar fractal, corresponds to a phenomenon with an unlimited capacity for spatial dispersion and with an a priori undefined variance.
In a fractal Brownian motion model (Eq.6), the Hurst exponent is allowed to vary from 0 to 1 (Huang and Bradford, 1992).A fBm is an expansion of the random walk or Brownian model (Bm) characterized by Hurst exponent H =0.5, first proposed by Mandelbrot and van Ness (1968).Conversely, the random walk or Bm model can be considered an especial case of the fBm.
For a fractal transect the semivariance function, according to Mandelbrot (1983, p. 353), is a function of the spatial separation: The fractal dimension, D, of a fractal surface or profile represented by its semivariogram can be estimated from the slope of the straight line portion of the semivariance, γ (h), versus the lag distance, h, when plotted on a double logarithmic scale.
For a two dimensional data set of point heights, the fractal dimension, D SMV , is computed from the Hurst exponent, H , obtained by Eq. ( 7), and the Euclidean dimension (d=3) as: In addition to the commonly used fractal dimension, D, parameter, the fractal roughness model requires a second parameter to define the relative position of the straight line of the variogram plotted on a log-log scale.Huang and Bradford (1990) defined this parameter as the crossover length, l, and it should be used together with D for characterizing soil surface roughness at small distances.The semivariance function of the fMb fractal model may be described as a function of both, crossover length, l, and the exponent, H , as: In Eq. ( 9) γ (h)=h 2 when h=l, which explains the origin of the term crossover length, first proposed by Huang and Bradford (1992) for soil microrelief quantification.The crossover length, l, may be estimated from the slope of the straight line portion of a variogram by: where a, is the intercept of the straight line of the semivariance log-log plot at the y-axis.It follows that characterizing soil surface roughness by a fractal fBm model requires two parameters, fractal dimension, D, and crossover length, l, as shown by theoretical considerations (Huang and Bradford, 1992) and by experimental results (Huang and Bradford, 1992;Eltz and Norton, 1997;Miranda 2000;Vidal et al., 2005).

Statistical and geostatistical roughness indices
From the various indicators that appear in the literature to characterize soil surface microrelief, the statistical index random roughness and the two indices based on geostatistical concepts, limiting difference and limiting slope, were selected.Random roughness index was chosen because it is the most widely used microrelief index and also because of its performance for modelling water storage in soil microrelief depressions (Kamphorst et al., 2000).The reasons to select LD and LS were: first that LD index is based on elevation differences at large spatial intervals, and second, that LS index can provide insight into the configuration properties of the soil surface at small scale intervals.Tables 1 and 2 list RR, LD and LS values estimated throughout the growth period of maize and oats, respectively.Both RR and LD were sensitive to vertical differences in roughness between experimental plots.Taking into account the four tillage treatments, RR varied between 3.75 mm and 14.49 mm during the maize season and ranged from 4.56 mm to 14.97 mm during the oats season.Limiting difference values tended to be slightly higher than RR values as they fluctuate from 4.24 mm to 18.11 mm and from 5.34 mm to 17.04 mm in the maize and oats growth periods, respectively.
The initial RR values for surfaces with 0 mm rain in the maize season were 14.49, 14.25, 13.00 and 5.36 mm in BS, CT , MT and N T treatments, respectively.A similar trend was observed for LD initial values.Therefore, initial surface roughness evaluated by indices RR and LD in the maize season were similar for BS, CT , and MT treatments.Hence, under maize, the NT treatment had the lowest initial RR and LD values (Table 1).In the oats season, however, RR values were 14.97, 8.02, 7.40 and 12.12 mm for BS, CT , MT and N T treatments, respectively.So, under oats, the BS treatment had higher initial values of the vertical component of roughness as estimated by RR and LD indices than CT and MT treatments (Table 2).The effect of tillage on initial surface roughness depends on the tillage tool, the number of passes, and on other minor factors such as tractor speed, clay content and soil water content.In our case study, additional microrelief levelling during manual sowing of oats may explain the divergences between initial roughness of BS and CT treatments.Also, amount and management of previous crop residue may influence initial soil roughness.Noteworthy, in the N T treatment under oats values of RR, LD and l indices were high throughout the crop period.Moreover, with N T , the initial roughness under oats was more than two-fold higher than in the respective treatment under maize.This illustrates the effect of large amounts of previous crop residues, which may contribute to the formation of remarkable roughness on the soil surface.
The relationship between RR and LD is shown in Fig. 1.In our study the correlation coefficients between these indices were 0.944 and 0.943 for the 20 data sets acquired during the maize season and the 16 data sets acquired in the oats season, respectively.Both relationships were linear and significant (P <0.01).Linden and van Doren (1986) and Bertuzzi et al. (1990) also found a linear relationship between these two indices.
The limiting difference, LD=1/a, in Eq. ( 4a) is the asymptote value of the first-order semivariance, i.e. the sill of the first-order semivariogram.Indeed, RR corresponds to the square root of the sill of the second order semivariogram.Consequently strong correlation coefficients between RR and LD are expected.However, distinct regression equations in the maize and oats seasons in our study together with previous results (Linden and van Doren, 1986;Bertuzzi et al., 1990) suggest that there is not a general relationship between these two indices, in spite of its strong dependence for a given specific experimental condition.Although RR and LD are indicators of the height distribution of surface microrelief, they do not account for the spatial component, i.e. mutual location of higher and lower points.Spatial structure of the microtopography is critical for a thorough characterization of the configuration properties of the soil surface and for depressional storage evaluation.
Linden and van Doren (1986) stated that LS is soil surface slope at small intervals, because Z/ h would approach LS when h approaches zero.Therefore, LS should give information on the side slopes of structural units, such as large aggregates or clods, at small intervals.Moreover, on an idealized soil surface, maximum roughness depends on the side slope of the structural elements protruding the soil surface.The magnitude of LS was small when compared with those of RR or LD, as this index ranged from 0.013 to 0.264 in the oats season and from 0.022 to 0.142 in the maize season.In our study cases no significant correlation was found between LD or RR and LS.
Because of the linear nature of Eq. ( 3), LD and LS are essentially independent of the minimum sampling interval, which represent an advantage for analyzing data sets with various sampling grids.However, maximum soil surface roughness as assessed by the LD or the RR indices is independent of the size of the structural elements at the soil surface, which means that two different surface configurations may result in the same LD or RR values (Merril et al., 2001).This continues to be a major problem in characterizing soil surface microtopography.

Fractal parameters: fractal dimension and crossover length
The main results of fractal analysis during maize and oats growth periods are listed in Tables 3 and 4, respectively.These include: sampling date, cumulative rainfall, fractal dimension and its standard error, crossover length and its standard error, coefficients of determination, upper cutoff of the straight line portion of the semivariogram plotted on a loglog scale and number of data points in this first segment of the semivariogram.
The first portion of the structural function γ (h), when plotted on a log-log diagram, showed a similar trend in the 36 data sets studied, indicating the existence of a correlation between semivariance and scale at small scale intervals.An example is shown in Fig. 2. The graphed results show a straight-line portion of the semivariogram at short lag distances with a step slope and then a second portion, which could be approximately modelled by a horizontal sill, so that in this segment correlation between the structural function and distance in general is absent.
A self-affine model may quantify the first straight-line portion of the semivariogram.Thus, stable estimates of fractal parameters, D and l, could be obtained only from the first segment of the structural functions, γ (h), before the scale breaks in slope.This break in scale is mainly related with the size of the structural units, aggregates or clods, at the soil surface, consistent with previous work on soil surfaces recorded by pinmeter (Miranda, 2000;Vidal Vázquez et al., 2005, 2006) and by laser scanning (Huang and Bradford, 1992;Eltz and Norton, 1997;Davis and Hall, 1999;Vidal Vázquez et al., 2005).
Fitting the first straight-line portion of the structural function, γ (h), has been recognized as a critical step in fractal analysis when a self-affine fractional Brownian model is  used.In our case study the lower cutoff limits of self-affinity were about the same magnitude as the distance between the pinmeter needles, i.e. 30 and 20 mm for data sets acquired during maize and oats seasons, respectively.The upper cutoff limits varied between 82.6 mm and 215.2 mm in the maize season and between 82.6 mm and 151.8 mm in the oats season (Tables 3 and 4).The upper limits in this study were similar to those of a previous case study, where the semivariogram algorithm was applied (Vidal Vázquez et al., 2007).However, the upper cutoff limits were lower than those previously estimated with the root-mean-square (RMS) algorithm (Vidal Vázquez et al., 2006, 2007).
Because the distance at which scales break, when the γ (h) structural function is used, approximately matches the characteristic size of the larger clods, residue fragments or in general structural elements on the soil surface, this distance has also been regarded as a fractal parameter of considerable interest.This scale has been referred to as the correlation length (Vidal Vázquez et al., 2006).
In general, linear relations between structural function and scale covering at least two orders of magnitude are required for estimating the fractal dimension, D, with low standard error values (Miranda, 2000).But in our study cases, as before quoted, the range of self-affinity was between a lower cutoff of about 20 to 30 mm and an upper cutoff of about 82.6 to 215.2 mm.A millimetric grid resolution acquired by laser scanning would enhance the straight line portion of the semivariogram on a log-log plot against the scale (Miranda, 2000;Vidal Vázquez et al., 2005).Furthermore, in each of the 36 surfaces studied, the ratio (l 2 / l 1 ) between the upper (l 2 ) and lower (l 1 ) cutoffs of the structural function, γ (h), largely exceeds 2 1/D , which is the minimal condition to accept an experimental D value over a range of fractal selfaffinity (Pfeifer and Obert, 1989).
Accuracy in fitting a power law will depend essentially on the number of data couples of the straight line portion of the structural function, γ (h), versus scale, h.The results listed in Tables 3 and 4 show that the linear regression was fitted with 3 to 7 and 3 to 5 couples of data in the maize and oats seasons, respectively.Vidal Vázquez et al. ( 2006) used the root-mean-square algorithm (RMS) to estimate fractal parameters, D and l, from a sampling grid with 225 elevation points and 784 cm 2 in surface and the straight line portions of the structural functions were also described by a comparable small number of data couples.However, a 1.82 m 2 sampling grid with 3025 height data points yielded at least 7 data couples to estimate the straight line section of two structural functions, root mean square, RMS, and semivariance, γ (h) (Vidal Vázquez et al., 2007).Therefore, the main limitation of fractal analysis for assessing changes in soil surface roughness may be the number of elevation data points that are recorded in each observation.Large data sets are required for a good approximation in fitting the power law from which fractal dimension, D, and crossover length, l, are estimated.Indeed, the available couple of data for fitting the power law in the first straight line part of the γ (h) structural function increased when data sets measured by non-contact laser scanning with a resolution in the order of millimeters were available (Miranda, 2000;Vidal Vázquez, 2005).Standard errors of fractal dimension and crossover length calculated by the semivariogram method are also listed in Tables 3 and 4 for data sets acquired during the maize and oats seasons, respectively.Standard errors in estimating D varied between 0.02 and 0.25 under maize and between 0.02 and 0.26 under oats.Standard errors in estimating l ranged from 0.26 mm to 4.21 mm and from 0.10 mm to 10.44 mm under maize and oats, respectively.Therefore, standard errors in crossover length may be as high as its estimated values, or even higher.Vidal Vázquez et al. (2006) analyzed data sets recorded by pinmeter with a comparable small size, i.e. 225 height readings per plot, using the RMS algorithm, and found D standard errors being in the range from 0.008 to 0.023 and l standard errors being in the range from 0.006 mm to 1.040 mm.These results indicate that the RMS algorithm reduces errors in fractal parameters estimation when compared with the semivariogram algorithm by a factor of the order of one magnitude.Consequently, the use of square sampling grids instead of rectangular ones is recommended for data sets with a small number of data points, as those recorded by pinmeter.
Coefficients of determination for the straight-line portion of the γ (h) structural function were between 0.816 and 0.996 and between 0.542 and 0.978 for data sets acquired during the maize and oats seasons (Tables 3 and 4), respectively.The results obtained with the RMS algorithm, from 225 height readings, by Vidal Vázquez et al. (2006) were more precise, as the respective coefficient of determination varied between 0.972 and 1.000.
Fractal dimension values ranged from 2.53 to 2.87 and from 2.55 to 2.98 in the maize and oats seasons, respectively.Therefore, the 36-microrelief data sets studied showed antipersistent features (D>2.5), also in accordance with previous results of fractal parameters estimated in random microrelief data sets, obtained after correction for slope and tillage marks (Miranda, 2000;Vidal Vázquez et al., 2005, 2006, 2007).
The crossover length values estimated by the γ (h) structural function varied from 0.69 mm to 9.49 mm and from 0.44 mm to 11.59 mm under maize and oats, respectively.The magnitude of l values is also consistent with previous works on data sets acquired under field conditions (Vidal Vázquez et al., 2005, 2006, 2007).These results clearly indicate a larger variation in scale of the crossover length, when compared with the fractal dimension as maximum differences between experimental data sets of this later fractal parameter were of 0.31 units and 0.40 units under maize and oats, respectively.Therefore, values of crossover length show a much greater sensitivity to changes in microrelief than the fractal dimension.This reinforces the relevance of the crossover length parameter as a discriminator of vertical differences in roughness.
For a microrelief model based on fractal concepts, the significance of the crossover length should be emphasized as it allows differentiation between various degrees of soil roughness, whereas the most known fractal dimension would be an indicator of the spatial configuration of soil microtopography  ( Huang and Bradford, 1992).In quantifying soil surface roughness, the fractal dimension, D, can be taken as a relative measure of the spatial distribution of different size structural elements on the soil surface (Huang, 1998).However fractal dimension, D, does not provide information on roughness vertical component.Therefore, the fractal dimension as a descriptor of horizontal variations of soil roughness should be used together with an additional index describing differences in roughness height for a thorough evaluation of soil microtopography (Huang, 1998;Vidal Vázquez et al., 2006).The fractal parameter crossover length, l, and the geostatistical index limiting difference, LD, were compared, because both indices are thought to stand for the vertical component of soil surface roughness (Eltz and Norton, 1997;Huang, 1998).Figure 3 shows the relationship between l and LD.Correlation coefficients were 0.478 and 0.686 for the 20 data sets acquired in the maize season and the 16 data sets acquired in the oats season, respectively (P <0.05).
Again, LD defines the vertical component of soil roughness based on mean absolute elevation differences at relatively large distances by the sill of the first-order variance.However, l represents the intersect of the γ (h) structural function with the y-axis on a log-log scale.Therefore, the crossover length parameter is rather a measure of nugget effect or discontinuity at small distances, differing from the sill or variance, which gives vertical fluctuation statistics.In fact, the magnitude of the discontinuity at small distances depends on the vertical statistics, but it may depend also on the horizontal variation of soil roughness.The relationship between l and LD in Fig. 3 was distinct in the maize and oats seasons, pointing to differences in surface configuration between these two experimental periods.
Crossover length and fractal dimension values estimated during the two successive growth periods showed a significant correlation (P <0.05), as shown in Fig. 4. Therefore, the D value showed a trend, to increase as l increased, which is an expected result, given the dependence between these  parameters in Eq. ( 9).A similar relationship between D and l was also found in previous studies (Vidal Vázquez et al., 2006, 2007).Moreover, during the maize season l and D were significantly correlated (P <0.05) for each of the four study treatments as the respective correlation coefficients ranged from R 2 =0.71 to R 2 =0.81.Likewise, l and D were significantly correlated (P <0.05) in three out of four treatments during the oats season, as the respective correlation coefficients ranged from R 2 =0.71 to R 2 =0.90; in this case the exception was the CT treatment.So, different D and l values were associated with different soil tillage systems.Consequently, the couple fractal dimension and crossover length appear to be a pertinent descriptor of soil roughness.

Tillage, crop cover and rainfall effects on soil roughness decay
Overall, roughness indices RR and LD decreased during the growth periods of maize and oats, as a function of cumulative rainfall, in the tilled treatments, either left fallow, BS, or under vegetative cover, CT and MT , as shown in Tables 1  and 2. Crossover length in general also decreased with cumulative rain, as shown in Tables 3 and 4.However, the roughness decay during maize and oats seasons was faster in the BS treatment than in the CT and MT treatments.Furthermore, RR, LD and l indices, in general, were not substantially affected by cumulative rain in the N T treatment under maize and oats, whose surface was protected by previous crop residues.Therefore, because no changes are induced by rainfall in no-till systems with various quantities of crop residues, one single observation along the crop season allowed characterization of soil surface roughness.Regression equations were developed to evaluate the relationship between roughness decay as quantified by RR, LD and l indices and cumulative rainfall for BS, CT and MT treatments, presented in Figs. 5, 6 and 7.The independent variable was the ratio between the values of an index for a given cumulative rainfall amount relative to the value of that index for the initial surface, instead of the estimated row values listed in Tables 1 to 4. These equations had a quadratic or an exponential negative shape depending on the index and in all cases were fitted to honor the initial value for 0 mm rain.Therefore, the fitted general equation for quadratic and exponential roughness decrease were y=ax 2 − bx+1 and y=exp.(−ax),respectively, where y is the proportion of roughness relative to its initial value, x is cumulative rainfall and a and b are regression parameters.In the bare soil treatment, BS, left fallow, the RR index decreased 41% and 39% from initial values after 229 and 350 mm cumulative rainfall in the maize and oats seasons, respectively.The respective figures for the CT , MT and N T treatments under crop cover during the maize season were 60%, 51% and 87% and during the oats season they were 57%, 72% and 96%.Therefore, the roughness destruction evaluated by the RR index can be ranked as follows: BS>CT ∼ =MT >N T .(l), as a function of cumulative rainfall, under maize and oats with minimum ).The RR index decay was best fitted by a quadratic negative relationship given by RR t /RR 0 =aP 2 −bP +1, where RR t is random roughness for a given rainfall amount, P , and RR 0 is the initial random roughness.In the BS treatment a single quadratic relationship was fitted to all the data collected in the maize and oats seasons since the relative decrease in roughness as a function of cumulative rainfall agreed for both data series.In MT and NT treatments two different relationships were fitted to maize and to oats data series.The fitted equations showed that under MT and NT tillage treatments, RR decay was faster in the maize than in the oats season.The LD index decreased 41% and 52% from initial values after 229 and 350 mm cumulative rainfall in the maize and oats seasons, respectively, in the BS treatment, which was left fallow.The respective changes for the CT , MT and N T treatments under crop cover during the maize season were 57%, 57% and 105 % and during the oats season they were 55%, 73% and 95%.Therefore, the roughness destruction evaluated by the LD index can be ranked as follows: BS>CT ∼ =N T >N T , similar to that of the RR index.
The LD index decay was also best fitted by a quadratic negative relationship.However, better correlations were obtained in general between cumulative precipitation and RR than between cumulative precipitation and LD.With the BS tillage treatment LD decay was also modelled by a single regression equation for the maize and oats data series.The decay of LD with the MT and NT tillage treatments, however, was modelled using two different curves.The LD decrease was also faster in the maize than in the oats season.
Crossover length showed in general a well-defined trend, to decrease in the tilled treatments during the maize and oats seasons.In most of the rain sequences the decline of l from its initial value (0 mm rain) at the end of the season was very strong.For example, in the bare soil treatment, BS, l index decreased 16% and 4% from initial values after 229 and 350 mm cumulative rainfall in the maize and oats seasons, respectively.However, in the maize season, MT treatment started with a relatively low value of l parameter (4.73 mm) for 0 mm rain at sowing time, which increased after 35 mm cumulative rainfall to about double the initial value (9.49mm).This is an inconsistent result and will be discussed below.In this particular case, the value of l for 35 mm rainfall was not taken into account for regression purposes.
In general, l values in the BS, CT and MT treatments exhibit a rapid decline from the initial reference state (0 mm rain) during the earlier rainfall events.The exception was again the MT treatment in the maize growth period.Therefore, l decay as a function of cumulative rainfall, P , was best fitted by a negative exponential relationship given by the equation l t /l 0 =exp(−aP ), where l t is the crossover length for a given cumulative rainfall amount and l 0 is the initial crossover length for 0 mm rain.Crossover length was sensitive to differences in soil roughness conditions, allowing a description of microrelief decay due to rainfall in the tilled treatments, although better correlations between cumulative rainfall and the indices RR and LD were most commonly obtained (Figs. 5,6 and 7).
In a review of tillage and rainfall effects on roughness decay, Zobeck and Onstad (1987) described Random Roughness degradation caused by rainfall with the equation RRp=RR 0 exp (-0.026P) and reported that this equation explained 76% of the variance among 418 data sets for different tillage operations and soils.Eltz and Norton (1997) conducted experiments with simulated rain for measuring soil surface microtopography decline using a laser scanning device under fallow and soybeans.These authors found that the decrease in roughness as measured by the l index was more rapid in the earliest degradation stages than RR decline, after which changed very slowly.Moreover, Eltz and Norton (1987) described the RR and l decline versus kinetic energy of rain by quadratic and exponential relationships, respectively.
The values of the LS parameter in general were much lower at the end of the experimental period, i.e. after 229 mm rain in the maize season and 350 mm rain in the oats season, than at the initial soil surfaces with 0 mm rain.This notwithstanding, in most of the rain series LS values first showed a tendency to increase and then decrease with cumulative rain.An increase in the LS index indicates increased slope of structural units in the soil surface.Such slope increase may be the result of consolidation effects during the first rain events (Eltz and Norton, 1997).The LS decrease at the end of the growth period indicates reduction of slope at small distance intervals due to filling of small depressions around largest structural units and crust development.
There was little variation of fractal dimension with increasing cumulative rainfall.Small changes in D values versus cumulative precipitation were characterized by various patterns, so that no general trend was recognized.Values for D were similar to those obtained in previous studies by Eltz and Norton (1997), Huang (1998) and Vidal Vázquez (2005, 2006, 2007).It is noteworthy that mean values of fractal dimension in the maize season (2.707) were significantly lower (P <0.05) than those in the oats season (2.804), and this in spite of the fact that essentially the same tillage operations had been conducted in both periods.These differences in fractal dimension, thus, on configuration of the soil surface, may be attributed to factors such as contrasting soil water content and aggregate stability during tillage (Kamphost et al., 2000) or to plant growth effects (Martínez-Turanzas et al., 1997) which further modify the soil surface configuration.

Physical interpretation
On agricultural soils, traditional tillage by mouldboard ploughing, followed by chisel ploughing, creates the largest degree of roughness (Kamphorst et al., 2000).Other less dominant factors influencing the configuration of soil surface microrelief may be the number of passes of the tillage tool, i.e. primary or secondary tillage, soil water content and aggregate stability.
On bare soil, the main factors that cause roughness decay or levelling should be first, the kinetic energy of rain drops, and second, slaking associated to air entrapment effects or differential swelling after sudden wetting.Moreover, soil roughness can either decrease or increase during rainfall, depending on both the surface condition and processes occurring on that surface.Surface breakdown processes due to kinetic energy and/or slaking are likely to reduce soil roughness, whereas erosion processes lean towards roughness increase because of rill formation (Huang, 1998;Darboux et al., 2005).Since these processes may occur simultaneously, with one or the other dominating at different spatial locations, the net result mainly affects the rate of change of surface roughness during rainfall.In our case study, no erosion symptoms were observed on any of the tilled treatments, not www.nonlin-processes-geophys.net/15/575/2008/ Nonlin.Processes Geophys., 15,[575][576][577][578][579][580][581][582][583][584][585][586][587][588][589][590]2008 even in the bare soil, BS, treatment, left fallow.Therefore, no secondary roughness increases by erosion will be taken into account.Soil cover by plant canopy and/or cover crop and no tillage prevents soil surface from raindrop impact, reducing the rate of roughness changes.However, plant growth may increase soil surface roughness (Martínez-Turanzas et al., 1997) while modifying the configuration of the soil surface microrelief, mainly due to interactions between soil surface and root system.
Again on bare soils, and when slaking is neglected, the decay or degree or destruction of soil roughness should depend on the initial roughness and on the kinetic energy of raindrops.Therefore, under these conditions it may be assumed that roughness decay will follow first order kinetics, which means that the diminution of roughness per unit rain amount (or per unit kinetic energy of rain) depends on the degree of roughness that is still available.As a first approach, this can be expressed as: where R t , is the roughness at a certain time, t, R 0 is the initial roughness and k is a constant describing the soil susceptibility against roughness destruction.Our results mainly support the above physical interpretation.Roughness decreased with increased rain at different rates in the BS, CT and NT treatments.However, roughness decay did not always follow negative exponential kinetics.Various reasons may explain the experimental behavior of roughness decay in our case studies, which were fitted by exponential and quadratic functions, depending on the roughness index used (Figs. 5,6 and 7).First, if slaking occurs, it should result in a faster roughness breakdown than can be expected from the kinetic energy alone.Second, crop and residue cover modifies the rate of decay in roughness.Also, less dominant factors such as soil water content during tillage may play a role in roughness decay dynamics.Furthermore, maize and oats root systems may modify the soil surface configuration, hence, the roughness decay dynamics in different ways.
In some instances, roughness indices may slightly increase with the first rainfall after roughness was increased by tillage.This initial small increase in roughness was detected for example in the CT treatment during the maize season, where the RR index increased from 14.25 to 14.40 mm for rainfalls of 0 and 35 mm, respectively (Table 1).Also, in the CT treatment during the oats season, RR and LD detected slight roughness increases between the initial stage with 0 mm rain and the stage after 50 mm cumulative rain (Table 2).These effects have been previously reported (Eltz and Norton, 1997;Huang, 1998).Consolidation of loosened soil particles within the largest voids without significant reduction in largest clods size may originate a denser soil surface with greater roughness than the freshly tilled soil surface and also may lead to increased side slopes, as mentioned above for the LS index.
However, the big increase in crossover length with the MT treatment in the maize season (Table 3) when the initial soil surface and the stage after 35 mm cumulative rainfall was compared was more difficult to interpret.For a cumulative rain of 35 mm, the ratio l/l 0 can be considered as an outlier, because of the relative high value of the numerator.Then, the value of l 0 in the freshly tilled soil surface may have been underestimated, which would change the shape of the regression between l/l 0 and cumulative rainfall.
A major problem of the characterization of roughness in our study was the small size of data sets acquired with low technology devices.The randomization process of locating successive small measurement plots for characterizing a given treatment may partly explain the dispersion of the regression functions developed for roughness decay.This notwithstanding, our microrelief data sets have been found to be self-affine in a small range of scales with an upper limit matching the size of the largest structural units in the soil surface and a lower limit equal to the horizontal resolution of point heights measurements.Therefore, accuracy and reproducibility of the roughness indices and fractal parameters could be increased in different ways: i) a large measuring grid could be used and a higher number of readings could be taken, ii) replicate microplots could be measured on each plot and on each date, and iii) the same microplot on each treatment could be used on each date.
The fractal parameters D and l have been useful to further quantifying the soil surface configuration and to discriminate between soil uses and crop cover.Mean D values were higher in the oats season, which means a more rugged soil surface, even if parameters accounting for vertical roughness, RR and LD, did not follow this trend.Also, the relationships between RR and LD versus l in the maize and the oats periods were not equivalent (see Fig. 3), which is indicative of differences in soil surface configuration properties between the two crop canopies.

Conclusions
Microrelief measurements taken on small plots and consisting of 225 height measurements showed spatial correlation after slope and tillage trend removal in a limited range of scales, allowing to estimate two fractal indices, fractal dimension and crossover length.
Roughness decay with increasing rainfall was highest for bare soil left fallow and it was not noticeable for notill soil protected by both plant canopy and crop residue cover.Roughness decay of conventional tillage and minimum tillage under vegetative cover was in-between those of bare soil and no-till.Vertical roughness decay showed a trend to be faster under maize than under oats.
Both random roughness and limiting difference decreased from initial values in a similar proportion with increasing rainfall.Crossover length index,exhibited a faster decrease from the initial values than random roughness and limiting difference.Even thought there was little variation in fractal dimension with cumulative rainfall this parameter was found to be significantly higher in the oats than in the maize season.

35Figure 2 .
Figure 2.An example of the relationship between the structural function K p a log-log diagram and the scale.

Fig. 2 .
Fig. 2.An example of the relationship between the structural function γ (h) plotted on a log-log diagram and the scale.
39parameter (l), as a function of cumulative rainfall, under maize and o conventional tillage (CT).

Fig. 8 .
Fig. 8. Roughness indices (RR and LD), and crossover length parameter (l), as a function of cumulative rainfall along the maize and oats season for no tillage (NT ).

Table 3 .
Sampling date, cumulative rainfall, fractal parameters (D and l) with respective standard errors (S.E.), coefficients of determination (r), upper limit of self-affine behaviour (U.l.) and number of data couples for fitting the first straight line portion of the semivariogram (n) during the maize season for different tillage treatments.(BS=bare soil, CT =conservation tillage, MT =minimum tillage, NT =no-tillage).

Table 4 .
Sampling date, cumulative rainfall, fractal parameters (D and l) with respective standard errors (S.E.), coefficients of determination (r), upper limit of self-affine behaviour (U.l.) and number of data couples for fitting the first straight line portion of the semivariogram (n) during the oats season for different tillage treatments.(BS=bare soil, CT =conservation tillage, MT =minimum tillage, NT =no-tillage).