Nonlinear Processes in Geophysics Assessing soil surface roughness decay during simulated rainfall by multifractal analysis

Understanding and describing the spatial characteristics of soil surface microrelief are required for modelling overland flow and erosion. We employed the multifractal approach to characterize topographical point elevation data sets acquired by high resolution laser scanning for assessing the effect of simulated rainfall on microrelief decay. Three soil surfaces with different initial states or composition and rather smooth were prepared on microplots and subjected to successive events of simulated rainfall. Soil roughness was measured on a 2 ×2 mm2 grid, initially, i.e. before rain, and after each simulated storm, yielding a total of thirteen data sets for three rainfall sequences. The vertical microrelief component as described by the statistical index random roughness (RR) exhibited minor changes under rainfall in two out of three study cases, which was due to the imposed wet initial state constraining aggregate breakdown. The effect of cumulative rainfall on microrelief decay was also assessed by multifractal analysis performed with the box-count algorithm. Generalized dimension, Dq , spectra allowed characterization of the spatial variation of soil surface microrelief measured at the microplot scale. These Dq spectra were also sensitive to temporal changes in soil surface microrelief, so that in all the three study rain sequences, the initial soil surface and the surfaces disturbed by successive storms displayed great differences in their degree of multifractality. Therefore, Multifractal parameters best discriminate between successive soil stages under a given rain sequence. Decline of RR and multifractal parameters showed little or no association. Correspondence to: E. Vidal Vázquez (evidal@udc.es)


Introduction
The surface of agricultural soils is built up from clods and aggregates arranged in a complex system of successive macro and/or microstructures.Soil surface roughness (SSR) is defined as the configuration of soil microrelief at small scales, of less than one meter (Allmaras et al., 1966;Linden and van Doren, 1986;Huang and Bradford, 1990;Merril et al., 2001).Soil surface roughness, taken on a scale ranging from cm to mm, plays a very important role in increasing water infiltration and the amount of crop water available and reducing runoff on cultivated lands (Podmore and Huggins, 1981;Armstrong, 1986;Kamphorst et al., 2000).Besides, SSR has been demonstrated to influence water infiltration, splash amount, overland flow and runoff routing (Govers et al., 2000;Römkens et al., 2001;Gómez and Nearing, 2005), to reduce runoff velocity and, thus, to decrease soil detachment and transport (Cogo et al., 1983) caused by water erosion.
In agricultural fields, tillage implements break-up soil and produce abrupt roughness increases.Subsequently, rainfall gradually decreases roughness.Römkens and Wang (1986) identified four different components of soil roughness.Type 1 is microrelief variations due to the size of individual grains or aggregates.Type 2 is due to cloudiness, as a result from tillage break-up.Type 3 consists of systematic differences in elevation produced by farm implements, such as furrows.Type 4 represents topographic variations at the plot scale and at higher scales, such as slope at field or basin scale.Roughness of types 1 and 2 are referred to as random roughness, whereas roughness of types 3 and 4 are called oriented roughness.Notice that rills and gullies may also create oriented roughness when concentrated erosion occurs.
Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
E. Vidal Vázquez et al.: Multifractal analysis of soil surface roughness Oriented roughness due to tillage tools exhibit periodic patterns and is easy to quantify by geometric models.However, random roughness assessment implies a quantification of the spatial distribution of microrelief elements with various that are randomly oriented, which has been more challenging (Huang 1998a).This notwithstanding random roughness (RR) is calculated as the standard error (Allmaras et al., 1966) or as the standard deviation (Kamphorst et al., 2000) of point elevations relative to a reference plane.Therefore, this index represents a statistical measure of vertical topographic variations, implicitly assuming that there is no spatial variation in surface roughness (Eltz and Norton, 1997;Kamphorst et al., 2000).
The relationship between rainfall energy and microtopographic features controls many transfer processes on and across the soil surface, including overland flow and erosion.The effects of soil roughness are mainly evident during the early stages of erosion processes, when soil detachment by raindrop impact, seal development and runoff generation are major factors and surface flow and concentrated flow are absent or limited (Mitchell and Jones, 1978;Helming et al., 1998).In these stages, depressional storage affects incipient runoff and soil detachment by raindrop impact due to the presence of pounded water.One important effect is the higher infiltration rate of the soil with greater roughness, although the influence of microrelief tends to decrease or even disappear due to surface sealing as rainfall progresses (Moore and Singer, 1990).Depressional storage also reduces runoff, and the importance of this effect also decreases with cumulative rainfall.Therefore, the usual observations reported in the literature and the common manner that models account for roughness effects on erosion indicate that a smooth surface generally yields more runoff and sediment than does a rough one (Zobeck and Onstad, 1987;Renard et al., 1997;Hairsine and Rose, 1992).On the other hand, soil erosion may be lessened under conditions of greater roughness not only from the reduction of runoff but also due to the greater level of hydraulic resistance that dissipates the flow energy (Einstein et al., 1951).Spatial configuration of soil microrelief is decisive for a description of depressional storage areas filled with water and for assessing the connectivity of runoff pathways (Onstad, 1984).
Erosion processes, i.e. sediment detachment, transport and deposition interact on both the micro-and macroscales (Foster, 1982).On the microscale feedback and interaction may occur with the structural units responsible for random and oriented roughness.On a soil surface with random roughness features, detachment on the peaks of the largest clods and aggregates produces a local sediment load that exceeds the local transport capacity causing deposition in small depressions.Detachment reduces clod height and deposition fills local hollows, decreasing both the vertical roughness component and the depressional storage through time (Helming et al., 1998;Kamphorst et al., 2000).As clods and aggregates are eroded and depressions are filled the gap between detachment and deposition decreases, so that deposition ends as the soil surface becomes smooth.Although detachment is not selective, entrainment and specially deposition are highly selective, resulting in coarse particles being left behind at the soil surface and fine particles being deposited in depressions, which may cause severe topography changes at small scales.Therefore, the total flow energy is unavailable to transport sediments (Foster, 1982;Abrahams and Parsons, 1991).In addition, the resistance of soil to detachment affects erosion by raindrop impact due to the modification of the clod size distribution (Römkens and Wang, 1986;Moldenhauer and Kemper, 1969).
All of the former effects tend to cause erosion on rough soil surfaces to be less than that on a corresponding smoother one, for a given slope.Nevertheless, these effects are mainly active during the early stages of rainfall, although they are also perceptible until runoff fills and interconnects the soil depressions, and a complete drainage network is developed.Most studies on the effect of the roughness induced by tillage have been focused on these early stages of rainfall and have been conducted on small laboratory or field plots, where a full drainage network can hardly build up.However, when the effect of the initial roughness on runoff and soil losses was studied at a scale large enough for the overland flow to reach a natural velocity, the experimental results have challenged the conventional view previously commented, showing little effect of roughness on runoff and overall higher sediment yield on the rougher surface at steep slopes (Helming et al., 1998;Gómez and Nearing, 2005).Moreover, these experiments showed that during the initial stages of rainfall, runoff and erosion were delayed on the rougher surfaced soils.But, as the experiment progressed, both runoff and erosion were less affected by the roughness treatments.At the end, the total runoff amounts did not vary as a function of the roughness treatment, and the total erosion rates were actually greater for the rougher soil surface treatments.This may be explained in terms of the surface features being of the same order of magnitude or larger in size than the flow depth, which affected the spatial distribution of the overland flow (Abrahams and Parsons, 1990) and induced a higher degree of flow concentration on the rough surfaces (Römkens et al., 2001).
On the other hand, the near-surface hydraulic gradient, that is, drainage and seepage, has been shown to significantly affect erosion (Bryan and Rockwell, 1998;Huang and Laflen, 1996;Huang, 1998b;Owoputi and Stolte, 2001).Darboux and Huang (2005) conducted a laboratory experiment to assess effects of soil surface depressions on runoff initiation, water runoff, and soil loss under different subsurface moisture regimes (seepage and drainage) and upstream flow conditions (with or without runoff).During the experiment, depressions delayed runoff initiation by storing water in puddles and enhancing infiltration.Once an apparent steady state was reached, surfaces with initial depressions slightly increased water flux compared with initially smooth surfaces.This effect occurred for both drainage and seepage conditions and persisted even after the surface storage capacity became low.Then, the authors concluded that the results showed that roughness had no significant effect on particle flux and concentration both under drainage and seepage conditions.
The above literature review indicates that the roughness effect on erosion can be further compounded by surface and subsurface factors and processes occurring at different scales.Surface microtopography influence surface processes like erosion, deposition, infiltration, etc. and many of these surface processes also causes a change in surface morphology (Huang, 1998a).It follows that microrelief is the product of several feedbacks and multiscale interactions involving a complex spatial and temporal variability.
Soil surface roughness has been the subject of an increasing number of studies since the pioneering work of Burwell et al. (1963) and Allmaras et al. (1966).Soil surface roughness has also been described using fractal models (Armstrong, 1986;Huang and Bradford, 1992;Gallant et al., 1994;Eltz and Norton, 1997;Huang, 1998;Davis and Hall, 1999;Vidal Vázquez et al., 2005, 2006, 2007).In general, for natural surfaces correlation distance has been found to be equal or smaller than the size of the greatest structural units, aggregates or clods on the soil surface, i.e. in the range from a few centimetres to a few decimetres (Huang and Bradford, 1992;Helming, 1993;Vidal Vázquez et al., 2005).Nevertheless, the effects of short-distance correlation on overland flow generation are thought to be negligible at the plot or microtopography (Darboux and Huang, 2005).Actually, most depression storage models assume a completely uncorrelated soil surface (Mitchell and Jones, 1978;Moore andLarson, 1979, 1978;Kamphorst et al., 2000).This may be the reason why, in soil microrelief studies, fractals have been relegated to a narrow range of scales and to specialized technical applications (Vidal Vázquez et al., 2005).
A fractal model refers to a set and can be characterized by a single parameter, such as the fractal dimension, D, while a multifractal model refers to a measure and can be characterized by a continuous spectrum of fractal dimensions usually referred to as the generalized dimension, D q .Multifractals are spatially intertwined fractals.Basic information on the multifractal concept and on procedures to characterize multifractals can be found in several books (Everstz and Mandelbrot, 1992;Falconer, 1997).Multifractal analysis has been recently used in many fields including soil sciences (Folorunso et al., 1994;Caniego et al., 2005;Bird et al., 2006;Dathe et al., 2006;Grau et al., 2006;Ibáñez et al., 2006;Roisin, 2007).However, very little information has been reported about the possible application of multifractal concepts for characterizing soil microrelief (García Moreno, 2006;García Moreno et al., 2008).
The aim of this study is to describe the characteristics of soil surface microrelief decay under simulated rainfall employing multifractal concepts and to compare the widely used statistical index RR with the multifractal quantification.

Soils
Two medium textured soils were selected from Mabegondo (Coruña province) and Pastoriza (Lugo province), both located in Northern Spain, on the basis of differences in clay and silt fractions, organic matter content and structural stability.Next, these soils will be referred to as MA for Mabegondo and LU for Pastoriza.The soils were Umbrisols (FAO) equivalent with Inceptisols (Umbrepts) according to the US.Soil Taxonomy (Table 1).Briefly, the Mabegondo soil had been continuously cropped to corn (Zea mays, L.) in summer and left fallow in winter for several decades, whereas the Pastoriza soil had been under corn and winter cereal (Lolium multiflorum, L.) rotation.Organic matter content of the later soil was much higher than that of the former.Silt and sand content differences were also remarkable, with higher silt values and lower sand values for the Mabegondo soil.Aggregate stability was significantly greater in the soil from Pastoriza than in the soil from Mabegondo.

Microrelief data sets
Roughness measurements were performed in laboratory conditions.Soil surfaces were prepared using air dry aggregates from the top layer of the studied soils and packing them in small containers or trays horizontally disposed.The largest aggregates were 20-40 mm in diameter.The depth of the soil in the containers was 0.05 m.A sand layer below the artificial soil layer allowed free drainage.Aggregates with the largest diameters were randomly located on the soil surface, avoiding sorting.The initial soil surface was gentle leveled before starting laser scanning.
These initial conditions were reconstructed to simulate a natural seedbed so that random roughness was rather low.Disturbed situations were obtained by simulated rain.Two different rainfall simulators were used.The Mabegondo soil (MA4 and MA6) was subjected to simulated rainfall under a drop-forming device at intensity of 30 mm h −1 .The Pastoriza soil (LU1) was subjected to simulated rainfall produced by a nozzle system at intensity of 65 mm h −1 .
The soil surface was scanned periodically, before starting simulated rain and after successive cumulative rainfall applications.Elevations were measured with an automated laser relief meter.Sample spacing, or distance between points along a transect and between transects was 2 mm and vertical resolution was 0.1 mm.The configuration of soil topography was 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.

Generalized fractal dimension
The concepts of fractals (Mandelbrot, 1983) and multifractals (Everstz and Mandelbrot, 1992;Falconer, 1997) and its application to soil science (Pachepsky et al., 2000) and more specifically microtopography (Eltz and Norton, 1997;Huang, 1998a;Vidal Vázquez et al., 2005, 2006;García Moreno, 2006;García Moreno et al., 2008) has been well described in the literature, so we shall not reiterate it here.In our work the scaling of point elevation measurements is directly assessed to obtain multifractal parameters.
If a profile of point heights measurements in a two dimensional space is covered by boxes of side length δ, the number of such boxes, n (δ) needed to cover the experimental transect when δ→0 varies as: The fractal dimension, D, can be obtained by counting the number n of boxes required to cover the object under investigation for increasing box sizes δ and fitting the slope of a log-log plot.
Using the box-counting technique for estimation of a single fractal dimension implies that each box employed for covering a transect is counted regardless of the proportion of the area occupied with pixels of a given height class.In other words, monofractal calculation does not account for the mass contained in each box, so that all of them have the same weight.
However, the generalized dimension calculated using the box-counting (BC) method essentially reflects the mass contained in each box.Let us consider that a domain of size R 2 as the support of the measure, δ.To assess the heterogeneity of the measure, δ, the unit initial square L×L is partitioned into n boxes of size δ×δ, by successive divisions in dyadic scaling down.Thus, the intervals for downscaling are logarithmically spaced.
Multifractal methods can resolve highly complex patterns of arrangement of the point elevation measurements defining soil surface microrelief.This complexity is characterized by using a mesh of square boxes of side length δ so that to each region of the space the corresponding quantity µ(δ) can be assigned.In practice, to implement the multifractal analysis of a distribution supported on a plane, a set of different meshes with cells or subintervals with equal length is required.The measure µ i (δ) distributed over an interval of sizes was computed from the experimental point elevation data, h i .First H i are normalized, H i =(h i / i h i ), i=1, 2, ..., n, with n i=1 H i =1, and then the measure µ i (δ) assigned to those blocks was calculated by adding all contributions h i inside a box.
In our case the study regions were square areas of initial sizes 460 mm for MA data sets and 550 mm for LU data sets and the size of each cell was 2 mm.Thus, depending on the data set the range of scales in the available data varied from 1 to 230 cells and from 1 to 275 cells, respectively.
For each box ith the probability distribution is: where α i is the Hölder exponent characterizing density in the ith box.For multifractal measures, the number N δ (α) of boxes of size δ here the probability has values in the interval obeys a power law: In practice, using the box counting method, for every box, i, the probability of "containing object", also called the partition function, is obtained for different moments q which can vary from -∞ to +∞.For multifractal distributed measures, the partition function scales with the block size as follows: A log-log plot of the quantity χ (q,δ) versus δ for different values of q yields: where τ (q) is the mass scaling function of order q.Also τ q can be written as: Multifractal measures are primarily characterized by their spectrum of dimensions.The concept of generalized dimension, D q , corresponds to the scaling exponent for the qth moment of the measure.Based on the work of Rényi (1955) generalized dimensions are defined as: Therefore, τ and D q are related as: For a monofractal, D q is a constant function of q, so no additional information is obtained by examining higher moments.However, for multifractal measures, the relationship between D q and q is not constant.In this case, the most frequently used generalized dimensions are D 0 for q=0, D 1 for q=1 and D 2 for q=2 termed, respectively, capacity, information (Shannon entropy) and correlation dimension.
The capacity or box-counting dimension, D 0 , is independent of the quantity of mass in each box; it is the scaling exponent of the number of non-empty boxes and takes into account the fact that the boxes are occupied or not.
The information dimension, D 1 , gives the probability of occupation of the ith box of size δ, without taking into account the way in which the measure is distributed within each box.Thus, D 1 provides a physical characterization indicating how heterogeneity changes across a certain range of scales and it is also related to Shannon entropy index.D 0 and D 1 take the same value if all the boxes have equal probability.Note that using Eq. ( 7) D 1 becomes indeterminate because the value of denominator is zero.Therefore, for the particular case that q=1, the following equation is used: The correlation dimension, D 2 , describes the uniformity of the measure values among intervals.The generalized dimension, D q , is more useful for the comprehensive study of multifractals.Differences between D q allow comparison of the complexity between measured soil microrelief data sets.
Commonly the degree of multifractality is assessed from the curvature of functions involved in multifractal analysis, i.e. singularity spectra f (α), mass exponent function τ (q) or generalized dimension, D q (e.g.Cheng, 1999).In this work the main properties of the multifractality will be described by a few parameters taken from the functions D q and τ (q).A sigma shaped D q spectra is taken as an indication that the measure is multifractal, whereas quasi-linear spectra are close to monofractals.Thus, in homogeneous structures D q are close, whereas in a monofractal they are equal.Notice that, for a monofractal distribution, values of D 0 , D 1 and D 2 become similar.Therefore, if a distribution has a tendency to multifractality it will be observed that D 0 > D 1 >D 2 .On the other hand, also the properties of the function τ (q), specially the local properties of τ (q) around q=1, have been found to be useful for describing multifractality (Cheng, 1997(Cheng, , 1999)).Following Cheng's work from the mass exponent function τ (q) the main properties of the multifractality of the measure can be characterized by the parameter τ (0)-2τ +τ , named here the multifractality index (MI).If MI<0 the measure corresponds to a multifractal, whereas MI=0 indicates a single fractal or a non-fractal.Other general indices derived from τ (q) as τ (q)−2τ (0)+τ (−q) are usually proportional to the MI (Cheng, 1999).

Evolution of the vertical roughness component
The evolution of random roughness (RR) with cumulative rainfall is shown in Fig. 1.This statistical parameter varied from 3.39 to 4.09 mm in the MA4 sequence, from 3.00 to 2.13 in the MA6 sequence and from 4.72 to 5.10 mm in the LU1 sequence.In natural conditions RR may vary approximately between 1 and 40 mm (Kamphorst et al., 2000).Consequently, the studied soil surfaces were rather smooth, as it is expected for seedbeds.The decay of roughness as a function of cumulative rainfall exhibits noticeable differences between the soil surface initially dry (MA6) and the two surfaces initially wet (MA4 and LU1).In the MA6 sequence random roughness after 50 mm cumulative rain was 71% of the initial value.However, in the MA4 and LU1 sequences roughness decay as a function of rain was negligible and final values were even somewhat higher than initial ones.The soil susceptibility again roughness breakdown depends mainly on kinetic energy of rain and on air entrapment or differential swelling  after sudden wetting.In the wet surfaces MA4 and LU1 air entrapment was avoided and therefore the main driving force of microrelief decay was raindrop kinetic energy.Moreover, during the experiment with simulated rainfall the soil surfaces were partly covered by a water layer acting as mulch that protected the artificial seedbed from the drop impact.Therefore, in our study cases, slaking by air entrapment caused a faster roughness breakdown of the initially dry surface MA6.Notice also that the small increase in random roughness with increasing rain may be the result of surface consolidation effects and aggregate rearrangements induced by rainfall, which have been reported before (Eltz and Norton, 1997).These effects have been also observed in MA6 surface when comparing the initial dry stage and the stage after 5.0 mm rain.3.2 Multifractal parameters of soil surface microrelief

Partition function
Values of the partition function χ (q,δ) have been estimated for the whole available box size range in steps of 2 k , 0<k<7.The log-log plots of the normalized measures χ (q,δ) versus measurement scales, δ, calculated with Eq. ( 5) were examined to find out whether the spatial pattern of soil surface microrelief obeys power low scaling.Figure 2 shows two selected plots for the MA6 surface, those of the dry initial stage and the subsequent stage after 5 mm rain.Figure 3 shows two more plots that correspond to the LU1 surface at the wet initial state and after cumulative 195 mm rain, respectively.For q>0 the partition functions of all investigated data sets showed a positive slope with a distinct linear behaviour.In general, for q≤0 the partition functions exhibited a negative slope and also a clear linear behaviour.However, in some cases, a deviation from linearity was observed at δ values close to unity and for moments q equal or close to -10.Visually, the most noticeable departure from the straight-line model was detected for q=-10 at δ values close to unity on the initial stage of the MA6 surface, as illustrated in Fig. 2. In other words, the largest grid length deviates most from the straight line, although there is also some curvature at the opposite end, i.e. the smallest grid length.Estimations of D 0 using 8 regression points resulted in values higher than 2.000, with no physical meaning.Therefore, the point corresponding at the largest grid length in Figs. 2 and 3 was excluded from the regression analysis.Thus, all the calculations were done choosing 7 regression points, so that the last point on the left was discarded.
One of the most important steps in multifractal analysis is to determine the range of δ and q exhibiting linear behaviour.Particularly for q<0, D q values may vary depending on whether all the regression points or only the points of the straight line region are used in the analysis.This issue is recurrent ever since multifractal analysis was first applied and a careful study of coefficients of determination is required (e.g.Evertsz and Mandelbrot, 1992;Bird et al., 2006;Grau et al., 2006, among others).Coefficients of determination, R 2 , of the straight line log χ(q,δ) versus log δ, standard errors of the slope, together with the corresponding D values for selected q moments, are listed in Table 2.All the calculations were done choosing 7 regression points as explained above.For q=-10, values of R 2 were higher than 0.999.For q=10, values of R 2 were higher than 0.992.It follows that for all the studied microrelief conditions and statistical moments (q=-10 to 10) the logarithm of the normalized measures versus the logarithm of the measurement scales fit a straight line with R 2 >0.992.
The distribution of a measure is considered as a monoor multifractal when the moments obey power laws, i.e. the double log plots of χ(q,δ) against log δ varies linearly.There were, however, differences in the degree of power law scaling between and within the three artificial surfaces submitted to simulated rain sequences.Moreover, for a given data set the linear fittings for moments q 1 and for q 1 often showed divergences.For instance, the coefficients of determination of the generalized dimension values, D max , calculated for q=10 were greater in all the surfaces of the LU1 sequence, than in those of the MA4 and MA6 sequences.However D min values calculated for q=-10 showed lower coefficients of determination in the surfaces of the LU1 sequence.Consequently, analysis of the partition function indicates by now different degrees of multifractality, which will be next addressed by analyzing the generalized dimension function, D q , and some derived parameters.
Results for D 0 , D 1 , and D 2 are listed in Table 2. Values for the capacity dimension, D 0 , were always 2.0 for the two soils, Mabegondo (MA) and Pastoriza (LU), indicating that the support of the measure is the Euclidean plane.The information dimension, D 1 , ranged from 1.998 to 1.984 and from 1.998 to 1.981 in the MA and the LU soils, respectively.The correlation dimension, D 2 , oscillated between 1.996 and 1.968 in the MA soil and between 1.997 and 1.961 in the LU soil.Coefficients of determination for D 0 , D 1 , and D 2 were equal to 1.000 both, for the Mabegondo and Pastoriza soils.Standard errors given in Table 2 are the standard errors of the slope obtained with linear regression.The errors for D 0 were equal to 0.0 in both, Mabegondo (MA) and Pastoriza (LU) soils.In the MA surfaces the ± deviation ranged from 0.000 to 0.003 and from 0.001 to 0.005 for D 1 and D 2 , respectively.In the LU surface the ± deviation for D 1 and D 2 oscillated from 0.000 to 0.002 and from 0.000 to 0.005, respectively.
Statistical moment, q, acts as a scanning tool scrutinizing the denser and rarer regions of the measure.For q 1, regions with a high degree of concentration are amplified, while regions with a small degree of concentration are magnified for q 1.Values of D for q=-10 and q=10, i.e.D min and D max , respectively also are shown in Table 2. Coefficients of determination for D min and D max were lower than those for the three first moments, as quoted above.For D −10 the highest standard error was 0.099, which corresponded to the MA6, surface with cumulative 55 mm of rain.Likewise for D 10 the highest standard error was 0.022 in the LU1 initial soil surface.
The value of the information dimension, D 1 , has been also considered as a good index of the heterogeneity in spatial distribution of a measure.The closer the D 1 value to the capacity dimension, D 0 , the more homogeneous is the distribution of the measure.In general, the width of the multifractal spectrum could be a practical parameter for characterizing and comparing soil surface roughness in microplots.However, the width parameter may be assessed in different ways.Frequently, either it has been referred to as the deviation of the D (q>0) from the D 0 values, as given by the ratio D min /D 0 , or it has been considered as the amplitude of the maximum and minimum dimensions, D(q min )−D(q max ), and, in this case, positive or negative q=10 are commonly retained.Notice also that uncertainties in estimating Dq values for q<0 lead to errors in the width parameter when estimated both by the D min /D 0 ratio or by the amplitude of positive and negative q max =10, i.e. (D min −D max ) as discussed by Tarquis et al. (2003).Therefore, the degree of multifractality depends also on the uncertainty, i.e. ± errors of these parameters.Standard errors are additive, and from Table 3 it follows that for (D −10 −D 10 ) they range between from 0.004 to 0.115.
In general for assessing multifractal parameters it should be taken into account that the higher the absolute q value Nonlin.Processes Geophys., 15,[457][458][459][460][461][462][463][464][465][466][467][468]2008 www.nonlin-processes-geophys.net/15/457/2008/ the worse the linearity of plots χ(q, δ) versus δ, and also the higher the uncertainty as measured by the standard errors.Consequently, in retrospective, the uncertainty of the width parameter when calculated by the amplitude of positive and negative q max =10, i.e. (D −10 −D +10 ), would be greater than the (D 0 −D 2 ) counterpart had been used.However, when extreme values are of interest indices defined on the basis of higher order statistical moments, q, may be advantageous.
On the other hand it has been reported (Cheng, 1997(Cheng, , 1999) ) that the MI index, earlier defined as τ -2τ +τ characterizes the main property of the multifractality of a measure.MI usually gives the best index in terms of minimum errors, because of the uncertainty increase as moment of order q increases.Notice also that from Eq. ( 8) it follows that τ (0)=-D 0 and τ =D 2 , but τ <>D 1 .For conservative multifractals, as in the case of the measure defined in our work, the MI becomes MI=-(D 0 −D 2 ).
As indicators of soil surface heterogeneity, and multifractal behavior Table 3 lists three roughness indices: (D min −D max ), D min /D 0 and the MI index.First, both, (D −10 −D 10 ) and D min /D 0 are positive for all the studied data sets.Second, the MI has negative values in accordance with the above discussion and with the convex property of τ (q) function at q=1.
The degree of multifractality in the LU1 sequence clearly decreased with increasing cumulative rain, as indicated by the four multifractal indices analysed (Table 3).In the MA4 sequence, with the soil surface initially wet, (D 0 −D 2 ), (D −10 −D 10 ), and D min /D 0 displayed and decreased trend during the first rainfall events and then further increased with cumulative rain, whereas (D 0 -2D 1 +D 2 ) approached zero.In the MA6 sequence with the soil surface initially dry, parameters, (D 0 −D 2 ), (D −10 −D 10 ), and D min /D 0 showed no a definite trend by increased rainfall whereas MI was <0 for the two last rainfall events.
A larger width of the D q spectra is associated with a higher heterogeneity of the soil microrelief features, whereas a decreasing trend in D q width could be regarded as a measure of homogenization.Similarly MI values closer to cero are indicative of a relatively low degree of multifractality.Initial soil surfaces are constituted by aggregates with a range of sizes relatively heterogeneous.All the three initial soil surfaces are characterised by quite large values of the parameters (D −10 −D 10 ) and D min /D 0. Likewise, these three initial soil surfaces exhibit MI values that are somewhat distant from cero.During the first rain events the transport capacity at the microscale is limited and the dominant processes causing soil surface disturbance are crusting and deposition.Crusting leads to vanishing of small-sized aggregates, whereas deposition of sediments produced by raindrop impact or eventually other mechanisms such as air entrapment, will reduce microrelief differences.Consequently, a more spatially homogeneous soil surface microrelief is observed.The lower degree of multifractality in the LU1 sequence with increased simulated rain matches these observations.How- ever, when erosion starts to overtake the deposition process, some of the previously deposited materials are removed and small micro-rilling may start, causing increasing heterogeneity at small scales.The trend to first decrease and then increase the degree of multifractality with increasing rain observed in the MA4 sequence as well as the rather opposite trend in the MA6 sequence may correspond to changes in the relative intensity of erosion and deposition processes by increased cumulative rainfall.

Vertical roughness decay and multifractal parameter evolution
The importance of taking into account, not only the commonly used fractal dimension parameter, D, but also a crossover length parameter, l, which gives insight into the vertical scale, when using the fractal approach, for characterizing differences in soil surface roughness between microplots, was emphasized by Huang (1998a).Besides, advantages of the joint use of these two parameters to quantify soil microtopography have been illustrated (Eltz and Norton, 1997;Vidal Vázquez et al., 2006, 2007).
Figure 5 shows the evolution of multifractal parameters (D min −D max ), D min /D 0 as well as the multifractality index, MI, defined by τ -2τ +τ as a function of random roughness, RR, for the three studied rain sequences.RR following successive rainstorms.However, the corresponding changes in width of the D spectrum, as assessed by (D min −D max ) or D min /D 0 , as well as multifractality index MI, allow to discriminate between data sets of a given rain sequence.In general, all of these three indices show a similar behaviour when plotted against RR, even if they are not proportional.
The LU1 sequence was the only where the values of (D min −D max ), D min /D 0 decreased whereas MI approached to cero as a function of cumulative rain (Table 3).In this case study, however, changes in RR were negligible.
(D min −D max ) and D min /D 0 varied from 0.632 to 0.032 and from 1.232 to 1.005 respectively, whereas MI changed between -0.0386 and -0.0032 when the initial (0 mm rain) and final (195 mm rain) soil surfaces were compared.These indices exhibits great differences between data sets measured at 0, 65 and 130 mm rain, but they were virtually equal at 130 and 195 mm rain.These results imply that, with increasing rain, the soil surface becomes more homogeneous as the parameters accounting for width of the D spectrum decrease and as the multifractality index comes close to cero.
The MA4 sequence also goes through minor changes in RR when submitted to successive storms.Parameters describing spatial heterogeneity, (D min −D max ) and D min /D 0 decrease and the MI varies from -0.0230 to -0.0043 with the first rain event of 10 mm at 30 mmh −1 intensity.The second event raises cumulative rain to 55 mm, but this results in rather slight changes in the value of the above parameters.However, increasing rainfall from 55 to 85 mm dramatically increases the heterogeneity of the soil surface, as shown by differences in the values of (D min −D max ), D min /D 0 and MI.This may be indicative of dominant erosion processes at the small scale during the last event of this rain sequence.
Finally, in the initially dry MA6 sequence, both RR and multifractal parameters change with increasing rainfall, but these fluctuations are unalike.The initial RR slightly increases after 5 mm rain due to consolidation at the soil surface and them steadily decreases wit cumulative rain.Afterwards values of (D 0 −D 2 ), (D min −D max ) and D min /D 0 exhibit a slight decrease by 5 mm rain and then they enlarge by 12.5 and 27.5 mm rain.Again, values of these parameters decrease by 50 mm rainfall.The MI parameter follows a comparable trend.This means that the spatial heterogeneity may be increasing or decreasing whereas the vertical component of roughness is decreasing.Imbalances, between erosion and deposition, dominant processes may be the cause of this type of soil surface evolution.
The above case studies indicate that there is hardly any correspondence between multifractal parameters, describing spatial heterogeneity or the degree of multifractality, and the statistical parameter random roughness, RR, which characterizes the vertical component of soil microrelief.These results are pointing out that the information registered at the horizontal and vertical scales by multifractal and statistical analysis, respectively, is complementary.For that reason, both RR and multifractal parameters should be taken into account when describing and modelling soil surface microrelief.

Conclusions
Multifractal formalism was appropriated for analyzing the variability of point heights measurements on a 2×2 mm 2 grid and could be a practical way of assessing the spatial heterogeneity of soil surface microrelief.Several multifractal parameters, such as (D min −D max ), D min /D 0 , and the multifractal index MI defined as τ -2τ may be useful for expressing the microrelief irregularities of the soil surface.The values of these indices were consistent with experimental observations and reflected different degrees of multifractality of the soil surface microrelief data sets.All of the three indices have been found to be sensitive to assess soil spatial heterogeneities, but the MI was the best index in term of errors in the calculations.
In general, changes in multifractal parameters obtained from the generalized dimension, D q for the mass exponent function, τ (q) under simulated rainfall showed no or little correspondence with the evolution of the statistical parameter random roughness, RR, which explains the vertical decay of soil microrelief.Multifractals were able to discriminate soil microrelief data sets with similar values of vertical component of roughness, thus accounting for the spatial configuration of microtopography.

Fig. 1 .
Fig. 1.Random roughness (RR) versus cumulative rainfall for three different initial soil surfaces, under given rain sequences.

Fig. 2 .
Fig.2.Plots on a log-log diagram of the partition function, χ(q,δ), versus measurement scale, δ of two soil surfaces from the MA6 sequence.

Fig. 3 .
Fig. 3. Plots on a log-log diagram of the partition function χ (q, δ), versus measurement scale, δ of two soil surfaces from the LU1 sequence.

Fig. 4 .
Fig. 4. Generalized dimensions and standard errors as a function of q for MA4, MA6 and LU1 sequences.

Fig. 5 .
Fig. 5. Width of the generalized dimension spectra as measured by (D min −D max ), D min /D 0 , and multifractality index, MI, as a function of random roughness, RR.

Table 1 .
Composition and main characteristics of the soil surfaces used in this study.

Table 2 .
Cumulative rain and multifractal parameters for three soil surfaces with different initial state conditions subjected to simulated rain.

Table 3 .
Cumulative rain, random roughness (RR), and roughness indexes derived from multifractal parameters for three soil surfaces with different initial state conditions subjected to simulated rain.These indices are: (D min − D max ), differences between generalized dimensions for q=-10 and q=+10, D min /D 0 , or ratio between D −10 and D 0 , and MI or multifractality index calculated as τ −2τ +τ .