A multifractal approach to characterize cumulative rainfall and tillage effects on soil surface micro-topography and to predict depression storage

Most of the indices currently employed for assessing soil surface micro-topography, such as random roughness (RR), are merely descriptors of its vertical component. Recently, multifractal analysis provided a new insight for describing the spatial configuration of soil surface roughness. The main objective of this study was to test the ability of multifractal parameters to assess in field conditions the decay of initial surface roughness induced by natural rainfall under different soil tillage systems. In addition, we evaluated the potential of the joint use of multifractal indices plus RR to improve predictions of water storage in depressions of the soil surface (MDS). Field experiments were performed on an Oxisol at Campinas, S ão Paulo State (Brazil). Six tillage treatments, namely, disc harrow, disc plough, chisel plough, disc harrow + disc level, disc plough + disc level and chisel plough + disc level were tested. In each treatment soil surface micro-topography was measured four times, with increasing amounts of natural rainfall, using a pin meter. The sampling scheme was a square grid with 25 × 25 mm point spacing and the plot size was 1350 × 1350 mm ( ≈1.8 m2), so that each data set consisted of 3025 individual elevation points. Duplicated measurements were taken per treatment and date, yielding a total of 48 experimental data sets. MDS was estimated from grid elevation data with a depression-filling algorithm. Multifractal analysis was performed for experimental data sets as well as for oriented and random surface conditions obtained from the former by removing slope and slope plus tillage marks, respectively. All the investigated microplots exhibited multifractal behaviour, irrespective of surface condition, but the degree of multifractality showed wide differences between them. Multifractal parameters provided Correspondence to: E. Vidal Vázquez (evidal@udc.es) valuable information for characterizing the spatial features of oil micro-topography as they were able to discriminate data sets with similar values for the vertical component of roughness. Conversely, both, rough and smooth soil surfaces, with high and low roughness values, respectively, can display similar levels of spectral complexity. Although in most of the studied cases trend removal produces increasing homogeneity in the spatial configuration of height readings, spectral complexity of individual data sets may increase or decrease, when slope or slope plus tillage tool marks are filtered. Increased cumulative rainfall had significant effects on various parameters from the generalized dimension, Dq, and singularity spectrum,f (α). Overall, micro-topography decay by rainfall was reflected on a shift of the singularity spectra,f (α) from the left side ( q 0) to the right side ( q 0) and also on a shift of the generalized dimension spectra from the right side ( q 0) to the left side ( q 0). The use of an exponential model of vertical roughness indices, RR, and multifractal parameters accounting for the spatial configuration such asD1 or D5 improved estimation of water stored in surface depressions.


Introduction
Soil surface micro-topography has been demonstrated to be a factor relevant for assessing the hydrological response of small plots (Darboux and Huang, 2005;Antoine et al., 2009) and also has been used to extract indicators of soil structure (Hairsine et al., 1992).In particular soil surface provides micro-catchments for rain, which affects runoff initiation during the first stage of the rainfall-runoff process.In soil science, surface roughness has been defined as a measure of variation in surface micro-topography at the square meter scale (Allmaras et al., 1966;Huang, 1998).
E. Vidal Vázquez et al.: Cumulative rainfall and tillage effects on soil surface micro-topography The predominant roughness index for studies of microtopography is the so-called random roughness (RR), which is calculated from the standard error (Allmaras et al., 1966) or from the standard deviation (Currence and Lovely, 1970) of point height readings from the mean value along a transect or on a surface.Prior to RR calculation, experimental height readings are corrected to adjust for the effects of slope and tillage tool marks.Various other indices exist for characterizing soil micro-topography (see Kamphorst et al., 2000 for a review).
Although there is not a standard filtering procedure before RR calculations are carried out, in general corrections of height readings involves two steps: (i) removal of slope effects and (ii) removal of tillage marks effects.Moreover, soil surface roughness is usually partitioned into oriented roughness and random roughness (Römkens and Wang, 1986).Oriented roughness is characteristic to a specific tillage tool and is obtained from experimental elevation data, after slope effects have been removed.Random roughness results from adjustment for both slope and tillage tool marks and consist of structural units of different sizes, e.g.clods, aggregates and grains, randomly oriented on the soil surface.Therefore, the "random" roughness concept has been developed to quantify features of structural elements with various dimensions on the soil surface, which are assumed to be randomly oriented.Notice that the use of the "random" term does not mean that height readings are randomly distributed, i.e. without spatial correlation (Huang, 1998).
Several devices have been developed to obtain point height readings, as recently summarized by Moreno et al. (2008) and García Moreno et al. (2008b).Briefly, until now the pin meter, which consists of a row of pins that can be lowered onto the study surface still remains the most commonly used equipment in field conditions.Indeed, the pin meter is a destructive method and, therefore, care should be taken to avoid possible disturbance of the surface; it also does not allow micro-relief measures over the same plot during successive periods with increasing cumulative rain.The relatively large grid spacing is another disadvantage of this method (Vidal Vázquez et al., 2005, 2006).Alternatively, non contact methods such as laser scanner (e.g.Huang, 1998;Vidal Vázquez et al., 2007), stereo-photography (e.g.Wagner, 1995), shadow analysis (García Moreno et al., 2008b), etc., allow to measure point elevation at higher resolution than pin meter and also eliminate surface disturbance.In spite of these disadvantages, pin meters are cheap, simple and reliable in field conditions.
On agricultural land, soil surface roughness is influenced by several factors such as tillage operations, vegetation, soil type and previous amount and intensity of rainfall (e.g.Allmaras et al., 1966;Zobeck and Onstad, 1987;Kamphorst et al., 2000;Vidal Vázquez et al., 2006;Paz-Ferreiro et al., 2008) and may be influenced by less dominant factors such as freeze-thaw cycles (Hansen et al., 1999).Soil tillage modifies surface roughness by breaking large clods into smaller ones and by introducing mounds, rips and furrows.Tillage operations result in abrupt changes in roughness, depending on its type and intensity.Rainfall produces decay of roughness left by tillage.Consequently, soil surface roughness is created by tillage and later smoothed by rainfall.
The water storage capacity of the soil surface following rain events, here referred to as Maximum Depression Storage (MDS), largely depends on soil micro-topography.Early studies recognized that rougher surfaces store more water than smoother surfaces and steeper slopes store less water than gentle slopes (Moore and Larson, 1979;Ullah and Dickinson, 1979;Onstad, 1984;Hansen et al., 1999;Govers et al., 2000;Gómez and Nearing, 2005).Thus, a smoothly tilled surface has little depression storage.More recently, it has been shown that the main effect of water stored in soil surface puddles was to retard the time to start runoff after rainfall begins (Darboux and Huang, 2005).Right after tillage a large roughness is often associated with a high infiltration capacity, which would result in a delay in runoff initiation at the beginning of rain events.This effect is more evident in soils where a surface seal develops.However, a higher surface roughness and, therefore, larger depression storage does not always reduce runoff and/or erosive soil losses (Helming et al., 1998;Darboux and Huang, 2005).
Some attempts have been made to measure directly MDS (Kamphorst and Duval, 2001), but most frequently it has been indirectly estimated from micro-relief measurements, assuming conditions of zero infiltration.Computer models based on geometrical algorithms have been developed for calculating MDS at the microplot scale from a grid of point elevation measurements (Moore and Larson, 1979;Onstad, 1984;Hansen et al., 1999).These algorithms identify individual depressions by determining local minima and then the stored volume is calculated by gradually filling such depressions to the overflow point.Also MDS estimations from Digital Elevation Models (DEMs) have been performed (Ullah and Dickinson, 1979;Huang and Bradford, 1990;Kamphorst et al., 2000) using depression-filling methods.
Depression storage has also frequently been obtained from empirical relationships with roughness indices (e.g.Onstad, 1984), the most often used being random roughness (RR), according to Kamphorst et al. (2000).In these studies, surface elevations were measured with different grid spacing and vertical resolution.On the other hand, roughness indices and depression storage have been often calculated by different methods from the topographic data sets.Therefore, results cannot always be compared.Notwithstanding, the main criticism to linear regression analysis between RR and MDS lies in that random roughness gives no information about the proportion of different aggregate sizes on the soil surface.If one parameter like RR is used to characterize the total surface of the plot, the spatial distribution of surface roughness will be disregarded.As a consequence, surfaces with the same RR may have a different MDS.
Several authors have recognized that mathematical description of the spatial characteristics of soil surface microrelief still remains a challenge, because most of the indices employed are merely descriptors of the vertical component of soil surface roughness (Kamphorst et al., 2000;Vidal Vázquez et al., 2005;Moreno et al., 2008).The widely used RR 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).Some attempts have been made in the past to account for the spatial component of soil surface roughness.Linden and van Doren (1986) developed the so-called limiting elevation difference (LD) and limiting slope (LS) to take into account the spatial aspect of roughness.LD and LS are based on the first order variogram of height differences.Huang and Bradford (1992) first proposed a fractal model to express soil roughness as a function of two fractal parameters, i.e. fractal dimension, D, and cross-over length, l.Fractal behaviour of soil surface roughness has been demonstrated over a limited range of scales (Huang and Bradford, 1992;Eltz and Norton, 1997;Vidal Vázquez et al., 2005).Crossover length best characterizes vertical component of surface roughness, while fractal dimension was considered as an indicator of the distribution of point elevations on the surface (Huang, 1998).
Multifractal theory permits the assessment of complex phenomena in a fully quantitative manner, for both spatial and temporal variation.Multifractal techniques and notions are increasingly widely recognized as the most appropriate framework to analyze the scale dependency of geophysical data, including topography, and also their extreme variability over a wide range of scales (Lovejoy and Shertzer, 2007).Although the application of multifractal theory to soil surface micro-relief is incipient, both, data sets acquired by pin meter (Moreno et al., 2008;García Moreno et al., 2008a) and laser scanner (Vidal Vázquez et al., 2008) have demonstrated multifractal behaviour.
In a previous work (Vidal Vázquez et al., 2008) we carried out a multifractal analysis of micro-topography data sets obtained in laboratory conditions by laser scanner under simulated rainfall.However, in many countries devices measuring point elevations at the mm scale are still not available.On the other hand, low technology devices have been frequently employed to measure microrelief with a resolution of the order of cm, and it has been proposed to take advantage of them for an effective evaluation of roughness parameters in field conditions (Merril et al., 2001;Vidal Vázquez et al., 2006).Thus, we perform a field experiment to measure micro-topography on a Brazilian Oxisol by pin meter.The choice of this low technology device was driven by practical reasons.Duplicate data sets were obtained for six different tillage treatments, first just after tillage and subsequently with increasing cumulative rainfall.As shown in our previous work (Vidal Vázquez et al., 2008) random roughness Table 1.General properties of the topsoil, at the 0-20 cm depth.pH H2O pH KCl O. M. (%) Sand (g kg −1 ) Silt (g kg −1 ) Clay (g kg −1 ) 5.9 5.5 3.79 410 80 510 evaluation merely allows description of the vertical component of soil micro-topography.Thus, firstly we focused our work on spatial characterization of the point elevation data sets obtained in field conditions.Moreover, in this work, we studied the experimental oriented and roughness random conditions, which was not taken into account before.Even thought the experimental grid was relatively coarse, it was found to be suitable to perform multifractal analysis.In a further step we employed multifractal parameters to assess the effect of tillage system plus cumulative rainfall on soil surface microrelief in field conditions.In addition, the experimental data sets were also used to estimate surface storage capacity by a depression-filling algorithm.Therefore, the main objective of this study was to test the ability of multifractal parameters to assess in field conditions decay of initial surface roughness induced by natural rainfall under different soil tillage systems.Our specific objectives were: (1) to assess the usefulness of the multifractal approach to characterize oriented and random roughness features and (2) to evaluate the potential of the joint use of multifractal indices plus RR to improve predictions of water storage in soil surface depressions.

Site, soil, tillage operations and rainfall record
The field measurements for the present study were obtained at the Centro Experimental de Campinas, of Instituto Agronômico (IAC), Campinas, SP, Brazil, located at 22 • 53 S and 47 • 04 W, at an average altitude of 60 m a.s.l.(above sea level).The soil was a Latossolo Vermelho Eutroférrico típico (Oliveira et al., 1989), according to the Brazilian Soil Classification System, equivalent to an Eutrudox (Soil Survey Staff, 2006).As shown in Table 1, the topsoil (0-20 cm depth) had a slight acid pH (pH H2O = 5.9 and pH KCl = 5.5), the organic carbon content was 3.79% and the texture was clayey (41% sand, 8% silt and 51% clay).
Long-term mean annual rainfall at the study site is approximately 1380 mm, from which about 1050 correspond to the rainy season lasting from October to March during austral summer.Mean annual temperature is 20.5 • C.
The surface roughness measurements were made between October and November 2000, at a site in which the slope was 5.10 m m −1 .A total of forty-eight samples were obtained corresponding to six tillage treatments over four dates (the initial stage just after tillage and three successive periods www.biogeosciences.net/7/2989/2010/Biogeosciences, 7, 2989-3004, 2010 Even though two data sets were obtained for each combination of tillage treatment and cumulative rainfall, they are not considered in this analysis as replications, but rather, as two independent duplicate measurements for each situation.Notice that cumulative rain by the third and fourth sampling dates varied between tillage treatments.This was because of the small length of the dry period intervals between two consecutive rain events following typical rain patterns of the austral summer in the study site.

Micro-topography data acquisition and processing
Point elevation readings were taken with a pin meter consisting of a row of 54 probes placed in a frame and spaced at 25 mm (Fig. 1).Pins were designed to slide up and down until the tip just impact the soil surface.The position of each pin along a row was photographically recorded and thereafter digitized (Wagner and Yiming, 1991;García Moreno et al., 2008a).The profiles were registered by means of photographs using a digital camera.Image analysis was applied to obtain point heights.The Profile Meter Program (Wagner and Yiming, 1991) developed by the USDA-ARS Wind Erosion Research Unit of Kansas State University was utilized for this purpose.
The sampling scheme was a square grid with 25 × 25 mm point spacing and the plot size was 1350 × 1350 mm (≈1.8 m 2 ), so that each data set consisted of 3025 individual elevation points (Fig. 1b).
Height readings, obtained on natural soil surfaces need correction to adjust for the effects of plot slope and tillage  marks.In our work oriented and random roughness, were obtained by removing trends due to (i) plot slope and (ii) slope + agricultural practices, respectively (Vidal Vázquez et al., 2005, 2006).In the oriented roughness condition adjustment for slope was made simply using the plane of best fit for each microplot of about 1.8 m 2 surface.The residual surface after slope effect removal included both random roughness features plus periodic effects due to tillage.The random roughness condition is considered to be associated to the random disposition of aggregates and clods on the soil surface (Allmaras, 1966;Huang, 1998;Vidal Vázquez et al., 2005).Simultaneous adjustments for both linear and periodic trend effects associated to slope and tillage tool marks were performed by using a procedure proposed by Currence and Lovely (1970) which removes such effect by row and column.The average height value for the i-th row and j -th column is obtained from the experimental data through the equation: where Z ij is the corrected height at the i-th row and j -th column after removal of both the slope and the periodic trends, z ij is the original height reading for the i-th row and the j -th column, z i is the average for the i-th row, z j is the average height for the j -th column, Z is the average of all points.Therefore, micro-topography of each soil surface was described by three different data sets, namely point elevations experimentally obtained, and those corrected to represent the oriented and random roughness conditions.In each case, the soil micro-topography was single 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).Thus, from each experimental DEM two more DEMs were obtained, representing the oriented and the random roughness condition.
In this study, random roughness (RR) was calculated as the standard deviation of height readings after correction for plot slope and tillage tool marks.In addition, the standard deviations (SD) of experimental data sets without any transformation of the recorded elevation and those of data sets adjusted for the effects of plot slope (which describes oriented roughness) were also calculated.Both RR and SDs of experimental surfaces and surfaces with oriented roughness were calculated as follows: where Z(x i ) = height reading at location x i and n is the number of height readings.

Multifractal analysis of micro-topography
The multifractal analysis of a probability distribution on a rectangular (square) region of the plane requires a set of different grids with rectangular (square) cells.A common choice for the grids is to consider dyadic downscaling (Evertsz and Mandelbrot, 1992;Kravchenko et al., 1999).This may be achieved by successive dyadic partitioning of each side of the rectangular (square) region support of the measure with a factor δ = 2 −k (k = 1.2.3. . .).At each size scale δ, a number N (δ) = 2 2k of cells are considered and their respective measure are found from elevation data.Elevation data Z(x) should be normalized (Moreno et al., 2008;García Moreno et al., 2008 a;Vidal Vázquez et al., 2008) in order to have a probability measure: Where x i , x j are mean values at the i/j -th cell of the experimental grid.The probability distribution or mass distribution was then analyzed for multifractality using boxcounting based moment method (Halsey et al., 1986;Evertsz and Mandelbrot, 1992).This method involves estimation of three functions: mass exponent function, τ (q), singularity strength (or local scaling index), α(q), and singularity spectrum, f (α).
Briefly, according to Chhabra et al. (1989) a partition function, χ (q, δ), was first determined as follows: where moment q is a real number ranging from −∞ to ∞.
For multifractally distributed measures, the partition function scales with the grid size as: Therefore, the mass exponent function, τ (q), for each q value is obtained by plotting log χ (q, δ) versus log δ.
If the probability function µ i (δ) in the neighbourhood of the grid scales with the grid size as µ(δ) ∼ δ −α i , then, the singularity strength or Lipschitz-Hölder exponent α i characterizes scaling in the i-th cell or spatial region.Parameter α i , which next will be referred to as singularity exponent for short, is a crowding index quantifying the degree of concentration of µ: the greater this value, the smaller the concentration of the measure will be and vice versa.It is in fact the logarithmic density of the i-th cell of the partition of characteristic size δ (Feder, 1988).
The number, N (δ), of cells of size δ where the probability distribution has singularity exponents between α, and α + dα, is related to the cell size as N α (δ) ∝ δ −f (α) (Halsey et al., 1986;Chhabra et al., 1989).The function f (α) is a scaling exponent of the cells with common α and can be considered as the generalized fractal dimension of he set of cells with singularities α.
According to Halsey et al. (1986) and Chhabra and Jensen (1989), the scaling exponent, f (α), can be determined from the mass exponent function, τ (q), by Legendre transformation as: www.biogeosciences.net/7/2989/2010/Biogeosciences, 7, 2989-3004, 2010 E. Vidal Vázquez et al.: Cumulative rainfall and tillage effects on soil surface micro-topography and α(q) = dτ (q) dq (6b) A plot of f (α) vs. α is called a multifractal spectrum.It is a downward function with a maximum at q = 0. Multifractal spectrum quantitatively characterizes variability of soil parameters with asymmetry to the right and left indicating domination of small and large values, respectively.The width of the multifractal spectrum (w = α max −α min ) indicates overall variability (Moreno et al., 2008) similar to the nugget effects in geostatistics.
For each data set we calculated multifractal spectrum with q values ranging from −10 to +10 in steps of 0.5.This was fine enough to show the multifractal behaviour in the studied moment range.
Multifractal measures can also be characterized by their spectrum of generalized dimension, D q , defined by Grassberger and Procaccia (1983), based on the work of Rényi (1955) for all D q = 1, by the following scaling relationship: where χ (q, δ) is the partition function.For the particular case where q = 1 Eq.( 7) becomes indeterminate, so D q is estimated by l'Hôpital's rule.The generalized or Rényi dimension, D q , is a monotonic decreasing function for all real qs within the interval [−∞, + ∞].When q < 0, χ emphasizes regions in the distribution with less concentration of a measure, whereas the opposite is true for q > 0 ( Chhabra and Jensen, 1989).The generalized dimensions, D q for q = 0, q = 1 and q = 2, are known as the capacity, the information (Shannon entropy) and correlation dimensions, respectively.An object is monofractal if D(q) is constant; the equality D 0 = D 1 = D 2 occurs, however, only if the fractal is statistically or exactly self-similar and homogeneous (Korvin, 1992).Otherwise, if the measure has a tendency to multifractal type of scaling the following relationship holds: Also, the generalized dimension, D q , (Eq.7) is related to mass exponent function, τ (q) (Eq.5) by the relationship: If τ (q) behaves non-linearly with respect to q, the measure, m, is said to be multifractal.Notice that mathematically the multifractal can be completely determined only by the entire fractal spectrum function.A few characteristic values of the function, however, can be used to describe the main properties of multifractals.The properties of functions τ (q) and D q , specially the values at q = 0, 1 and 2 are frequently used for describing multifractality (Cheng, 1999;Vidal Vázquez et al., 2007).

Estimation of surface depression storage
Water stored on puddles of the soil surface, referred to as maximum depression storage (MDS) was estimated using an algorithm to fill depression on the soil micro-relief, similar to the way described by Onstad (1984).The estimation was simplified by not considering any infiltration.The algorithm first identifies local minima and, then, fills depressions stepwise from the bottom to the top.The search for the overflow point or outlet is performed by exploring the neighbourhood of each minimum.The number of neighbours used in our study was eight.
Boundary conditions have been demonstrated to be important for MDS calculation.Because of the small dimensions of the experimental plot, if depressions identified at the plot border were allowed to free drain (open boundaries), MDS values would be underestimated.The effect of depressions originated by random roughness and located at the plot border can be taken into account a method proposed by Kamphorst et al. (2000).Following this method, depressions running into the plot border and contributing to water storage can be regarded as features that continue outside the plot boundary and its volume calculated assuming a symmetrical pattern of micro-catchment distribution.In this work two boundary conditions were considered at the plot border: (i) free drainage and (ii) impeded drainage out of the random depressions at the plot boundaries.

Statistical analysis
SPSS, version 15.0 was used for all statistical analyses.An ANCOVA (analysis of covariance) was performed using tillage type (disc harrow, disc plough or chisel) and tillage intensity (primary tillage or two successive operations) as factors and rain as a covariable.This choice was employed because of cumulative rain intensity differences between treatments at the third and fourth measurement period (Table 2).Variables were log-transformed when it was required to meet the requirements of ANOVA (i.e.normality and heteroscedasticity of errors).A total of 48 samples (i.e. 3 tillage types × 2 tillage intensities × 4 measurement periods under increased cumulative rain × 2 duplicate samples) were analyzed per surface condition (i.e.experimental, oriented and random).Separation of mean values of multifractal parameters was subjected to the Tukey test.Correlations were made by Spearman ranking.

General characteristics of the soil surfaces
Figure 2 shows contrasted stages of the soil surface for two treatments, i.e. primary tillage by disc ploughing and disc ploughing plus levelling.In both treatments a freshly tilled soil surface and a disturbed surface, produced after 232.8 mm cumulative rainfall can be compared.Initial soil surfaces are permeable, whereas in rain disturbed surfaces a sedimentary or depositional crust, less permeable, was developed.Over the crusted surfaces small aggregates can be observed.They result from erosion of larger structural units such as clods and still had not been incorporated into the crust.Differences in micro-topography between the initial and the crusted stages are also apparent.Mean values of the standard deviation, SD, of the experimental date sets was 33.95 mm, whereas extreme values ranged from 16.03 to 57.40 mm.Removal of slope trend significantly (P < 0.05) reduced these figures, as for oriented roughness the mean of SD for height readings was 19.35 mm, being minimum and maximum 7.18 and 19.35 mm.In turn random roughness, RR, was significantly (P < 0.05) lower than oriented roughness, with a mean of 11.78 mm, and extremes ranging from 3.76 to 23.41 mm.Kamphorst et al. (2000) studied a number of surfaces with RR values between 1 mm and 44 mm, corresponding to relatively smooth seedbeds and vey rough conditions obtained by mouldboard ploughing.Therefore, our data set was mostly characterized by medium roughness conditions.
Mean values of roughness for the experimental and oriented conditions ranged as follows: disc plough > chisel plough > disc harrow > disc plough + leveling > disc harrow + leveling > chisel plough + leveling.

Multifractal analysis of experimental, oriented an random roughness
Multifractal analysis was carried out for a total of 144 data sets, because for each surface three different roughness conditions, i.e. experimental, oriented and random, were taken into account.The double log plots of the normalized measures χ (q, δ), versus measurement scales, δ, (Eq.5) were examined to find out whether the point elevation data obeys power low scaling.These plots also identify the range of moments that need to be considered to study the scale variation of point elevation data sets.Figure 3 shows selected examples of the partition function for oriented roughness and random roughness conditions, which was estimated in steps of 2 k , 0 ≤ k < 5. Visually, some noticeable departure from the straight line model was observed for moments q in Fig. 3b.Notwithstanding, determination coefficients, R 2 , were higher than 0.98 for the statistical moments in the range from q = −10 to q = 10, and this in all the 144 studied data sets.In general, the double log plots showed higher deviation from linearity for the highest q moments after removal of slope and tillage mark effects (i.e. for random distribution of clods, aggregates and grains on the soil surface), whereas in experimental and oriented roughness conditions they were practically linear (R 2 > 0.99).Data sets with a much larger number of height readings (> 50 000), measured by laser-scanning, allowed www.biogeosciences.net/7/2989/2010/Biogeosciences, 7, 2989-3004, 2010 Table 3. Summary statistics of selected multifractal parameters from the generalized dimension, D q , spectrum for experimental, oriented, and random surface conditions.
calculation of partition functions from eight dyadic steps 2 k (0 ≤ k < 7) and seven regression points (Vidal Vázquez et al., 2008).All of the studied micro-relief data sets in our study showed good scaling trends, even for a smaller range of k (0 ≤ k < 5).Therefore, as quoted by Moreno et al. (2008) a simple pin meter can provide valuable information for multifractal assessment of soil surface roughness, even if a laser scanner will produce more detailed information, which allows micro-topography analysis over a large scale range.The distribution of a measure is considered fractal (monoor multifractal) when the moments obey power laws (Evertsz and Mandelbrot, 1992).The scaling properties observed through analysis of the partition function can be further tightened up by determining if there is simple (monofractal) or multiple (multifractal) scaling types.
Figure 4 shows some examples of generalized dimension or Rényi spectra, D q .They correspond to rainfall stages with 0.0 and 120.9 mm, and have been drawn for experimental, oriented and random roughness conditions.Calculations were performed in the range −10 ≤ q ≤ 10.First, Rényi spectra from duplicated measurements taken in the same date can show noticeable different patterns.Second, most of the spectra were more or less sigma-shaped curves, which displayed more curvature for negative values of q than for positive ones; however spectra approaching the straight line model, which implies a monofractal scaling nature, were also recognized.Note that for a surface roughness with a monofractal tendency D q spectra would be quasi linear.All the sigma-shaped D q functions passed through 2.00 at q = 0 and matched minimum and maximum values at q 1 and q 1, respectively.
Table 3 lists summary statistics for several multifractal parameters cropped from the Rényi spectra for experimental, oriented and random conditions.A monofractal also would be characterized by D 0 = D 1 = D 2 .In all of the 144 surfaces D 0 > D 1 > D 2 , indicating that the soil micro-topography had a tendency to multifractal scaling property.Nevertheless, differences (D 0 − D 1 ) ranged from 0.002 to 0.067 and (D 1 − D 2 ) oscillated from 0.04 to 0.118, which suggests again various degrees in the homogeneity/heterogeneity of the soil surface roughness.In general the width of the D q spectra, assessed by parameters such as (D −5 − D 0 ) and (D −5 − D 5 ) indicated different degrees of heterogeneity in the associated measures, and this for experimental, oriented and random roughness conditions.
Determination coefficients, R 2 , were highest for moments q = 0 and q = 1 and decreased for lower and higher |q| moments.For q = 10, R 2 achieved values higher than 0.996, 0.997 and 0.999 in the sets of surfaces with experimental, oriented and random roughness, respectively.For q = −10, R 2 values for the corresponding surface conditions were 0.982, 0.980 and 0.983 (data not shown).Standard errors of D q values increased with increasing |q| values and they were much higher for left (q < 0) than for the right (q > 0) branch of the Rényi spectra.For |q| = 10, D q errors were much higher than for |q| = 5 (Fig. 4).Therefore, next the moment ranges −5 < q < 5 will be retained for discussion purposes.4. Generalized dimension, D q , spectra (−5 < q < 5) of surfaces tilled with disc plow, for experimental, oriented and random roughness conditions at two surface stages, 0 and 120.9 mm rain.
(Bars are standard errors).
The effect of cumulative rain on Rényi spectra, for the three different roughness conditions studied, is first described in Fig. 5. Here, each D q spectrum corresponds to a measurement date and so it was obtained as the mean value for two replicated measurements on neighbouring plots.Overall, features observed in Fig. 5 indicate various degrees of heterogeneity in the associated measures, depending on the cumulative rain.However, the random roughness condition was characterized by a much smaller degree of heterogeneity, as its mean D q curves under successive rain amounts showed a smaller variation in width than the respective curves for experimental and oriented conditions.In other words, some individual data sets from the three different surface conditions (experimental, oriented, random) studied here can exhibit gentle slopes of the D q vs. q function, which indicate rather Fig. 5. Mean values generalized dimension, D q , spectra (−5 < q < 5) of surfaces tilled with disc plow as a function of cumulative rainfall, for experimental, oriented and random roughness conditions.
homogeneous measures.On average, however, homogeneity increased in surfaces with random roughness conditions.Width of the D q branches, and, therefore, generalized dimension values for the studied q moments, increased or decreased between two successive measurement dates, so that they had no bearing with cumulative rainfall.On the other hand RR values steadily decreased with increasing cumulative rainfall.Therefore, changes in multifractal parameter values induced by rainfall showed no or little correspondence with the evolution of the vertical microrelief component, described by RR.On artificial surfaces, under simulated rainfall, Vidal Vázquez et al. (2008), found similar results.
The capacity or box-counting dimension, D 0 , was not significantly different from 2.00, which corresponds to a Euclidean support in the 144 studied surfaces.The entropy dimension, D 1 , quantifies the degree of heterogeneity of the distribution, µ itself by measuring the way its Shannon entropy scales as the linear size of the mesh shrinks.Parameter D 1 also has been considered as the average of the logarithmic densities or concentrations of the multifractal distribution.Therefore, D 1 may be viewed as the expected value of the different concentrations when the distribution itself is taken into account, and it also determines the geometrical size of the set where the "main part" of the distribution concentrates.The values of D 1 ranged from 1.933 to 1.998, 1.950 to 1.996 and 1.973 to 1.998 for experimental, oriented and random conditions, respectively (Table 3).Mean D 1 values were significantly higher (P < 0.05) for random roughness compared to experimental and oriented roughness.In spite of this, some data sets from the three different roughness conditions studied here were characterized by D 1 values close to 2 and Dq functions with rather gentle slopes.Again, this pattern was more frequent for random roughness than for experimental or oriented roughness.
The correlation dimension, D 2 , computes the correlation of measures contained in intervals of size δ.The values of D 2 ranged from 1.882 to 1.996, 1.906 to 1.992 and 1.952 to 1.996 for experimental, oriented and random roughness conditions, respectively.Also mean D 2 values were significantly higher (P < 0.05) for surfaces with random roughness than for surfaces with oriented or experimental roughness.
Table 3 also lists summary statistics for various parameters describing the width of the Rényi dimension, D q .On average, width of the D q spectra was significantly narrower for surfaces with random conditions than for those with experimental or oriented conditions.Parameters extracted from the Rényi spectrum for q < 0, i.e.D −5 , (D −5 − D 0 ), etc. and parameters given its total width such as (D −5 − D 5 ), were significantly lower (P < 0.05) for the random roughness condition, whereas no significant differences were found for parameters obtained for q > 0, such as D −5 and (D 0 − D 5 ).This suggests that removal of slope plus tillage trends mainly affects the left part of the Rényi spectrum, where q < 0.
The singularity spectrum is a powerful tool to analyze similarity or difference between the scaling properties of the measures and it also permits to asses the local scaling properties of individual point elevation data sets.First, the symmetry/asymmetry of the singularity spectrum is and indicator of homogeneity/heterogeneity. The wider the spectrum is, (i.e., the largest the α max − α min value) the higher the heterogeneity in the scaling indices of the measure and vice versa.Second, the branch length of the f (α) spectrum gives insight about the abundance of the measure.Small f (α) values at the end of a long branch correspond to rare events, whereas the largest f (α) value is the capacity dimension that is obtained at q = 0.
Selected examples of singularity spectra, f (α) vs. α are shown in Figs. 6 and 7.Those depicted in Fig. 6 correspond to experimental, oriented and random roughness for the disc harrow treatment, and those in Fig. 7  versus two successive tillage operations.Table 4 lists summary statistics of parameters obtained from the singularity spectra for the selected range of q moments (−5 < q < 5) and for three roughness conditions, i.e. experimental, oriented and random.All singularity spectra in Figs. 6 and 7 are characterized by a concave down shape, but they present very different patterns regarding symmetry features.When comparing singularity spectra from duplicated data sets measured on the   same treatment and date, thus with similar values of random roughness, RR, great differences may be observed because of the natural variability between neighbouring microplots.Most of the f (α) vs. α spectra showed an asymmetrical curve with a wider and also frequently a longer right side, where q < 0. However, this was not a general rule, as the left branch of the singularity spectrum was wider than the right one in 7, 8 and 10 out of 48 cases, for experimental oriented and random conditions, respectively.The magnitude of changes around the maximum value of f (α) is a measure of the symmetry of the singularity spectrum.So differences (α max − α 0 and α 0 − α min ) indicate the deviation of the spectrum from its maximum value (q = 0) towards the right side (q < 0) and the left side (q > 0), respectively.On the other hand, also in most of the studied cases the right branch of the f (α) spectrum was longer than the left one, but in several cases the opposite was true.Small f (α) values correspond to rare events.Therefore, a longer right branch of the f (α) spectrum suggests depressions of micro-topography were very rare, whereas a longer left branch suggests peaks of height readings were less frequent.Altogether, asymmetry toward the left side indicates domination of large or presence of extremely large values in the spatial variability of height Table 4. Summary statistics of the width of the singularity spectrum, f (α) vs. α for experimental, oriented, and random surface conditions.(Calculations were made for the moment range: −5 < q < 5).
(q −5 , q 5 ) α max α min α max -α min α max -α 0 α 0 -α min readings, while asymmetry toward the right shows domination of small or presence of extremely small values of height readings.
In summary, the most frequent pattern of the f (α) spectrum for experimental and oriented conditions was an asymmetrical curve with a greater tendency toward the right side, where q < 0, which indicated dominance of small point elevation values.This information is very important, since the values of q > 0 are directly associated with the measure that, in our case, is in turn related with the height reading.Small clusters with peaks of height readings embedded within larger areas with small height readings caused this pattern of asymmetry.The random roughness condition was characterized by a much narrower and symmetrical f (α) spectrum (Fig. 6).So mean widths (α max − α min ) for the range of q moments −5 < q < 5 were significantly lower (P < 0.05) for experimental and oriented conditions than for random conditions.The smaller mean values of the width, (α max − α min ) of the singularity spectra and the more symmetrical shape of the f (α) vs. α spectra of the random condition indicate a more homogeneous spatial distribution of height readings when compared with experimental and oriented conditions.
Mean α max and (α max − α min ) values, calculated in −5 < q < 5 range, were also significantly lower (P < 0.05) for the random condition than for the experimental and oriented conditions.However, mean α min and (α max − α min ) values were not significantly different (P < 0.05) within the three random conditions studied.These results corroborate the differential pattern of homogeneity/heterogeneity between the random condition and the experimental or oriented condition made evident by the generalized dimension, D q , analysis (see Figs. 4 and 5).Therefore, it may be said that, on www.biogeosciences.net/7/2989/2010/Biogeosciences, 7, 2989-3004, 2010 average, the random roughness condition, resulting from simultaneous removal of slope and tillage tool marks, is associated with increased homogeneity of the residual height readings.This is an expected result, because random roughness is due to the spatial distribution of elements such as clods and aggregates randomly oriented on the soil surface (Huang, 1998).Nevertheless, it should be taken also into account that according with our results, homogeneity of individual data sets may increase or decrease when slope or slope plus tillage tool marks are filtered, even if in most of the cases trend removal produces increasing homogeneity in the spatial distribution of height readings.
The effect of tillage intensity, i.e. primary tillage versus two successive operations, on the singularity spectrum is illustrated in Fig. 7, where disc ploughing and disc ploughing plus leveling are compared for the oriented roughness condition.All the 16 concave down curves depicted in this figure are characterized by a long right side of the f (α) spectrum.Although the mean vertical component of the micro-topography, RR, was significantly higher in the primary tillage, no significant differences (P < 0.05) were found in the mean widths, (α max − α min ), of the f (α) spectrum.In spite of some minor differences in their dimensions and symmetry, similarities at both low and high f (α) values on treatments with two tillage intensities, i.e. disc plough and disc plough + leveling, suggest analogous or related micro-topographical features.Thus, multifractal analysis shows that rough and smooth soil surfaces with high and low RR values, respectively, can display similar levels of spectral complexity and heterogeneity.Moreno et al. (2008) and García Moreno et al. (2008a) performed multifractal analysis by pin meter for three different tillage treatments in three different experimental fields.Each data set consisted of 10 000 points.The texture was sandy clay loam for two of the soils and sandy loam for the third soil.The width of the singularity spectrum were calculated in the moment range −5 < q < 5, so that w a = (α −5 − α 5 ).For the experimental, i.e. real roughness condition, w a oscillated from 0.09 to 0.99, 0.07 to 0.292 and 0.010 to 0.395 under chisel, roller and mouldboard plough, respectively.The values of w a for the random roughness condition ranged from 0.003 to 0.121, 0.003 to 0.034 and 0.003 to 0.019 under chisel, roller and mouldboard plough, respectively.Trend removal of micro-topography data sets also decreased width of the f (α) spectrum, with the exception of the chisel treatment.Even if García Moreno et al. (2008b) used a gliding box algorithm to perform multifractal analysis, singularity spectra were much wider in our clayey soil (Table 4), with mean values, w a = 0.462, 0.405 and 0.239 for experimental, oriented and random data sets, respectively.Soil texture, antecedent rainfall and other factor such as soil structure and organic matter content may have had an effect on the complexity of surface micro-topography, which would explain the above differences in width of the singularity spectrum.

Effects of cumulative rain, tillage type and intensity on multifractal parameters
Results of the ANCOVA analysis are summarized in Table 5.
Cumulative rain had significant effects on several multifractal parameters.Increasing rainfall showed a trend to increase D −5 , whatever the roughness condition (P = 0.004 for experimental, P = 0.058 for oriented and P = 0.076 for random roughness conditions).The greater effect of cumulative rainfall was observed on the total width, i.e. (D −5 − D 5 ) of the generalized spectrum, D q , (P = 0.002 for experimental, P = 0.039 for oriented and P = 0.020 for random conditions).Moreover, increased cumulative rainfall caused decreasing D 5 (P = 0.007, P = 0.052 and P = 0.002 for the experimental, oriented and random conditions respectively), and, therefore, the width of the right branch, q > 0, i.e. (D 0 − D 5 ) of the D q spectrum also decreased (data not shown).In addition, the capacity dimension, D 1 , and the correlation dimension D 2 also showed a trend to decrease with decreasing cumulative rain, so that the greatest effects were observed on the experimental data sets and the smallest ones on the oriented data sets.
Likewise Table 5 shows that cumulative rain had also significant effects on parameters describing the width of the singularity spectrum, f (α).Increased rain showed a trend to increase α −5 values of the singularity spectra (P = 0.008, P = 0.103 and P = 0.055 for experimental, oriented and random roughness conditions, respectively) and therefore, the width of the right branch (q < 0), i.e., (α −5 − α 0 ).Parallel to results obtained for D q , increased rainfall also decreased both, α 5 significantly (P = 0.038, P = 0.027 and P = 0.001 for experimental oriented and random roughness respectively) and (α 0 − α 5 ) (data not shown).Because the increase in the width of the right side, (α max − α 0 ), of the singularity spectra was larger than the decrease in the width of the respective left side (α 0 − α min ), the total width of the f (α) spectrum, (α −5 − α 5 ), increased with increasing rain (P = 0.006, P = 0.037 and P = 0.004 for experimental, oriented and random roughness, respectively).
The above results show that the main effect of increased rainfall on the multifractal characteristics of soil surfaces consists on a shift of the singularity spectra, f (α) from the left side (q 0) to the right side (q 0) as well as a shift of the generalized dimension spectra from the right side (q 0) to the left side (q 0).In both cases, these changes in symmetry/asymmetry indicate that the small values of height readings become more frequent and the largest ones become rarer as cumulative rainfall increases.On the other hand the increase in total widths (α max − α min ) of the f (α) spectrum with increasing cumulative rain indicate that the heterogeneity and, therefore, the spectral complexity of the soil surface was higher in the crusted stages than in the initial ones, just after tillage was performed.Therefore, multifractal parameters give a good description of what can be observed in different stages of the soil surface micro-topography (see Fig. 2).No significant effects (P < 0.05) were detected for tillage type, tillage intensity or the interaction tillage type × tillage intensity for the random roughness conditions.It is remarkable also that the interaction tillage type × tillage intensity showed significant effect for several multifractal parameters of the experimental roughness condition.

Prediction of maximum depression storage
Values of water stored in surface depressions (MDS) estimated taking into account the border effect were in the range from 1.82 to 8.77 mm m −2 , with a mean value of 3.67 mm m −2 .Kamphorst et al. (2000)    varying from 1 to 13 mm m −2 , for a RR range between 1 and 46 mm.The mean MDS estimated for free drainage conditions was 1.17 mm m −2 , about three times smaller than for the previous condition.Allowing free drainage of depression at the border of the microplot seems not to be realistic, because it is like consider a big pit surrounding a small surface, which would produce a important bias in MDS estimations.
Results of regression analysis between MDS and random roughness, RR, are shown in Table 6.Coefficient of determination was higher when the effect of microplot border was taken into account (R 2 = 0.679) than by free drainage conditions (R 2 = 0.551).The regression equation was: MDS = 0.333 RR − 0.248.For impeded drainage conditions, i.e., border effect, Kamphorst et al. (2000) obtained MDS = 0.28 RR, with a higher correlation coefficient (R 2 = 0.85).The above results show very significant correlation coefficients (P < 0.001).However, the dispersion of the data around the regression line is high, and therefore, prediction of MDS from RR are not satisfactory, even if the RR explains 85% of the MDS variance (Kamphorst et al., 2000).The main criticism to this lineal model remains in that RR describes only the vertical component of soil surface roughness.
Based on results of the multifractal analysis we tested exponential models that take into account parameters such as D 1 , or D 5 , which are thought to contain information about the spatial component of surface roughness.The general form of models was: MDS = a RR Dq + b and MDS estimations were performed with free and impeded drainage conditions.As shown in Table 6, the joint use of RR and D q somewhat improved MDS prediction.When MDS was estimated taking into account the border effect, the models with D 1 and D 5 as an exponent over the base RR, explained 71.5% and 73% of the variance respectively.Improved determination coefficients also have been obtained by the exponential models when MDS was estimated assuming free drainage at the border; again, D 5 performed better than D 1 .Although exponential models combining the vertical component of roughness, RR, and a multifractal parameter accounting for spatial features of soil roughness are still not accurate enough for a reliable MDS prediction, they do improve prediction from lineal models.These results warrant more research.

Conclusions
A field experiment was conducted to assess effects of cumulative rain and initial tillage treatment on soil surface roughness.Height measurements were obtained by pin meter at ≈1.8 m square meter scale.Remarkably good multifractal scalings were obtained when analyzing measures constructed from experimental micro-topography data sets and also from data representing the oriented and random roughness conditions, which were retrieved from the former after trend removal.On average, after slope plus tillage tool marks removal the spectral complexity of point elevations decreased.However, the spectral complexity of individual data sets can increase or decrease as a result of trend removal.
Data sets measured by a low technology device such as the pin meter have proven to be useful to perform multifractal analysis and, in turn, multifractal parameters showed advantages for assessing tillage and rainfall effects on soil microrelief.And this was in spite that spatial configuration patterns of soil micro-topography from duplicate neighbour measurements taken at the same treatment and on the same date can show great differences, for the three roughness conditions analyzed.
Similar levels of spectral complexity can be displayed by rough soil surfaces resulting from primary tillage and smooth ones produced by two tillage passes, with high and low random roughness values, respectively.Therefore, random roughness and multifractal parameters describe distinct microrelief features and produce different statistical information.
Increasing cumulative rainfall decreased the value of parameters from the right (q > 0) branch of the generalized dimension spectrum, such as D −5 , (D 0 -D −5 ), D 1 and D 2 , and increased the value of parameters from the left (q < 0) branch, such as D 5 and (D 5 -D 0 ).Likewise, increasing rainfall, decreased the width of the left (q > 0) side, (α 0 − α 5 ), while the width of the right (q < 0) side, (α −5 − α 0 ), of the singularity spectrum.This behaviour was observed for all the three roughness conditions studied i.e. experimental, oriented and random.
Depression storage (MDS) estimated by a depressionfilling algorithm was three times higher when depressions at the border were taken into account.Although results are still not accurate enough for a reliable MDS prediction, exponential models combining the vertical component of roughness, RR, and a multifractal parameter accounting for spatial features of soil micro-topography were more accurate than lineal models based solely on RR.

Fig. 1 .
Fig. 1. (A) Scheme of the pin meter used for point height measurements and (B) sampling grid in cm.

Fig. 3 .
Fig. 3. Selected plots of the natural logarithms of the partition function versus measurement scales.(A) oriented roughness and (B) random roughness.

Fig. 7 .
Fig. 7. Singularity spectra for disc plow and disc plow plus leveling under oriented roughness conditions.

Table 2 .
Summary of tillage treatments and cumulative rainfall (mm) at the four micro-topography reading periods.
compare primary tillage

Table 5 .
Effects of cumulative rain, tillage type, tillage intensity, and the interaction tillage type x tillage intensity on some multifractal parameters, for experimental surfaces and surfaces with oriented and random roughness.

Table 6 .
Regression equations to predict MDS from random roughness (RR), and random roughness plus multifractal parameters (D 1 and D 5 ).