Mapping soil fractal dimension in agricultural fields with GPR

We documented that the mapping of the fractal dimension of the backscattered Ground Penetrating Radar traces (Fractal Dimension Mapping, FDM) accomplished over heterogeneous agricultural fields gives statistically sound combined information about the spatial distribution of Andosol’ dielectric permittivity, volumetric and gravimetric water content, bulk density, and mechanical resistance under seven different management systems. The roughness of the recorded traces was measured in terms of a single numberH , the Hurst exponent, which integrates the competitive effects of volumetric water content, pore topology and mechanical resistance in space and time. We showed the suitability to combine the GPR traces fractal analysis with routine geostatistics (kriging) in order to map the spatial variation of soil properties by nondestructive techniques and to quantify precisely the differences under contrasting tillage systems. Three experimental plots with zero tillage and 33, 66 and 100% of crop residues imprinted the highest roughness to GPR wiggle traces (mean HR/S=0.15), significantly different to Andosol under conventional tillage ( HR/S=0.47).


Introduction
Numerous studies have documented scale invariance of soil and other porous earth materials over a broad range of scales (Oleschko et al., 2000;Caniego et al., 2005;Tarquis et al., Correspondence to: K. Oleschko (olechko@servidor.unam.mx)2006; Meng et al., 2006;Di Domenico et al., 2006;Jawson and Niemann, 2007).Self-similarity, a most striking property of isotropic fractals, means that each piece of a shape is geometrically similar to the whole (Mandelbrot, 1983).Self-affinity is a related concept referring to physical objects, time-series or recorded traces which must be scaled differently along the principal coordinates in order to conserve their shape.When a remote object with fractal nearsurface properties is explored with electromagnetic waves, the returned signals also become fractal.Microwaves, such as those emitted by the GPR (Ground Penetrating Radar) antenna, backscattered from the soil are self-affine functions of time whose fractal dimension is close to the mass fractal dimension of the high-permittivity, moisture-filled pores in the soil.Our group's use of the radar (Oleschko et al., 2002(Oleschko et al., , 2003) ) to find the soil's fractal dimension has provoked great attention, and was viewed by Jaggard (2002) as an extension of Richardson's idea of measuring the jagged cost of Britain with yardsticks of progressively smaller size (Mandelbrot, 2002)."In radar" -wrote Jaggard -"the wavelength λ plays the role of the yardstick, while the backscattered intensity counts the number of λ-sized scatterers".The pore topology, as well as the roughness of pore/solid interfaces are all imprinted in the GPR record's wiggle traces (Oleschko et al., 2003), this information can be decoded from the recorded wave trains using Mandelbrot's (2002) Fractal Geometry Toolbox.In the present research we selected for this "decoding" two dimension-estimators, Rescaled-Range statistics and the wavelet technique (Benoit software, SCION Corp., 2003).The mathematical model relating the GPR record's dimension to the mass fractal dimension of the Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.K. Oleschko et al.: Mapping of moisture in scale invariant soil high-dielectric-permittivity points in the soil had been developed and calibrated by our group (Oleschko et al., 2002(Oleschko et al., , 2003) ) using a combination of laboratory experiments, computer simulation, and field tests (Oleschko et al., 2003).
The present paper has three objectives: 1. To document in situ the GPR's ability to extract the fractal dimensions of the dielectric permittivity-and mechanical resistance soil spatial patterns; 2. To explain how this fractal dimension can be extracted from the GPR record; 3. To show that a proposed new way to combine the GPR traces' fractal analysis with routine geostatistics (kriging) makes possible to map the spatial variation of certain soil properties in fast and non-destructive mode.
2 Mathematical and physical background

Wave scattering on fractals
In the present research, we selected a widely accepted soil model (Pachepsky et al., 2000;Oleschko et al., 2002): a mass-fractal distribution of solids and pores, where the highpermittivity clusters are associated with water-filled pores.This model obeys the empirically observed scale invariance of the pore space (Korvin, 1992;Oleschko et al., 2000) and is simple enough to be treated analytically.If the solid grains and the pores are D-dimensional mass fractals then the nonuniform internal structure of an R-sized fractal aggregate manifests itself in nontrivial mass (M) and density (ρ) scaling (Zosimov and Lyamshev, 1995): and where D is mass fractal dimension, a 0 is a characteristic length of a single grain or pore, and E is the dimension of the embedding Euclidian space.The number of solid grains as well as pores of characteristic size R also follows a power law: In case of isotropy, the density correlation (4) The intensity of monochromatic waves scattered on a mass fractal (in our case on fluctuations of the refractive index) is proportional to where → q is the wave vector.If the correlation function of the refractive index fluctuations is proportional to that of the material density ρ( → x ), it follows from Eqs. ( 4) and ( 5) (Zosimov and Lyamshev, 1995) that, The observation of such power laws in scattering experiments provides both a verification of the fractal nature for the studied structure and a convenient way to measure its fractal dimension D. The use of Eq. ( 6) has become a common procedure to determine the fractal dimension of aerogels, colloidal aggregates, and polymers from small angle X-ray, neutron, or optical scattering measurements (Sinha, 1989;Schaefer et al., 1990;Jaggard and Jaggard, 1998;Huanling et al., 2006;Melnichenko and Wignall, 2007).

Estimation of the fractal dimension
Several types of dimensions are required for the complete description of a fractal; each one has its special significance and measures a specific set attribute (Korvin, 1992).The values of different fractal dimensions, or even the same fractal dimension estimated by different techniques, do not necessarily coincide.In the present study, the fractal analysis of the selected radar traces was accomplished using two techniques of the Benoit software (SCION Corp., 1999): the Rescaled-Range analysis and the wavelet analysis.In our previous research (Oleschko et al., 2002(Oleschko et al., , 2003)), both techniques proved to be sufficiently robust for microwave roughness measurements, yielding two fractal dimensions -D R/S and D Wand two corresponding Hurst exponents -H R/S and H W .The Rescaled-Range method divides the data set to "windows" of equal length w.For each window length the range R(w) and scatter S(w) are computed and the avarage value R/S(w)= R(w) S(w) estimated.For self-affine processes R/S(w) ∝ w H , that is the Hurst exponent H is found from the slope of the straight line that fits log R/S(w) vs. log w.The Hurst exponent is related to the fractal dimension of the graph as H =2−D.When H is close to 0.5, the behavior of the data is random (a total absence of persistence or memory).If 0<H< 1 / 2 , this indicates a negative persistence (antipersistence), while for 1 / 2 <H<1 the persistence is positive and increases when H changes from 1 / 2 to 1. Persistence is often identified with positive, and antipersistence with negative correlations.Positive correlation means that, if a positive (or negative) value is found at some position, with a high Nonlin.Processes Geophys., 15,[711][712][713][714][715][716][717][718][719][720][721][722][723][724][725]2008 www.nonlin-processes-geophys.net/15/711/2008/ probability a value with similar sign would be found at the neighboring site (Benoit, SCION Corp., 2003).
The wavelet transform method of estimating D has the advantage that it can be applied for nonstationary traces.In the Benoit software the mother wavelet is a step function, and n different wavelet transforms of the analyzed trace are computed using a i =2 i , i=0, 1, . .., n−1 as scales.Denoting by S i the variance from zero of the i-th wavelet transform, and calculating the numbers (Benoit, SCION Corp., 2003).Note that traces with H near zero are very rough (high amplitude oscillation), antipersistent, and have a fractal dimension close to 2, while traces with H near 1 are smooth (with homogeneous amplitude distribution), persistent, of fractal dimension near 1.

Spatial interpolation of irregularly sampled data (kriging)
In geostatistics (Robertson, 1998;Webster and Oliver, 2001) the spatial variability of a random regionalized variable f (x, y) is characterized by its semivariogram γ .Suppose we observe f at two points P (x, y) and Q(x + ξ, y + η), then γ is defined as In case of isotropy we have γ (ξ, η)≡γ ( ) where = ξ 2 +η 2 .For P ≡Q we have γ (0)=0.Sometimes the semivariogram does not go down all the way to zero when →+0.In this case the positive value ε 2 = lim →+0 γ ( ) is called nugget.If the regionalized variable f (x,y) is of zero mean value, homogeneous and isotropic, and the points P and Q are far from each other, Z(P )&Z(Q) become uncorrelated as . The value of this asymptotic flat region of the semivariogram is called sill, the -value, beyond which the semivariogram becomes flat is called range.Practically, the range is that distance within which the measured f -values are correlated, and beyond which, they become uncorrelated.The numerically found raw semivariogram is approximated by one of the following simple mathematical models (where σ 2 is the sill, a is the range, the separation distance, called lag, and ε 2 is the nugget): Spherical; Exponential; Gaussian; Pure nugget and Power Law models.Modern texts (Korvin, 1992;Hardy and Beier, 1994) use power law semivariogram models to describe fractal HC reservoirs.In this case this model has no sill, and no range.A special case of Eq. ( 7) is the linear model corresponding to α=1.
The main use of semivariograms is found in kriging (named after D. G. Krige, South African Mining Engineer), this is a procedure to find an unbiased optimal estimate for the unknown value of a spatial random variable X at point P if we know X P 1 , X P 2 , . .., X P N at N other points.For simplicity take N = 3, suppose we know X P 1 , X P 2 , X P 3 and we want to estimate X at a new point P as a weighted sum X P ,est ≈W 1 X 1 +W 2 X 2 +W 3 X 3 where the coefficients W i are selected in such a way as to minimize the expected errorvariance: V =E X P ,est −X P 2 = min subject to the condition W 1 +W 2 +W 3 =1.
The weights satisfy the kriging equation where γ (i, j )=γ P i −P j , is the appropriate model semivariogram at lag = P i −P j , λ is the Lagrange parameter.
3 Experiment description and data analysis

General considerations
To confirm experimentally the strong relation that exists between the mass fractal dimension of the pore space and the roughness of the scattered microwaves, we chose a volcanic soil (Andosol) with 1. well-documented self-similarity of both the solid and pore networks (Oleschko et al., 2000); 2. high microstructure stability, water retention capacity and infiltration rates (Oleschko and Chapa, 1989); and 3. slowly developed macrostructure, but extremely high biological activity (especially of moles, migrated to the plots with zero tillage) and exceptional erosion potential worsened by the geomorphology (Fig. 1).
We assumed that the pore topology and network continuity could be visualized on the georadargrams and are encoded in the roughness of the recorded wiggle traces.Our working hypothesis had been that if the georadargrams of two Andosol experimental plots with different moisture contents, resulting from contrasting management practices, are compared, the differences in the subsurface structure would affect the amplitude distribution of the scattered microwaves, and this could be decoded, by fractal analysis, from the recorded traces.The differences in soil subsurface structure are imprinted also in the roughness of GPR output frequency distribution (Fig. 2).

Site characterization
The studied site forms part of a long term experimental area started by CENAPROS (National Research Center for Sustainable Production of the National Research Institute of Forest, Agriculture and Animal Production -INIFAP) in 1995.This area is located southwest to the Patzcuaro Lake, Michoacan State, México (Fig. 1), and comprises two experimental fields each of 750 m 2 (50 by 15 m), and seven plots of 91.6 m 2 ( 4 by 22.9 m).All of them have a slope close to 9% (Fig. 1).In the last 12 years, detailed qualitative and quantitative analyses of the physical, chemical and biological properties of mollic Andosol under conservation or zero tillage (0, 33, 66 and 100% crop residues) have been accomplished, and the dynamics of crop yield was documented on the experimental fields under zero tillage and conventional management.Detailed dynamics of the physical, chemical, and biological properties of Andosol from the studied area under different management systems and crop rotations were described previously by Oleschko andChapa (1989), andTiscareño-López et al. (1999).Zero tillage (ZT) consists of minimum soil removal with a "zero tillage" machine seeding directly above the last year's residues, which covered about 38% of soil surface.
Conventional tillage (CT) comprises deep furrowing to a depth of 30 cm with a reversible 3-disk attachment disking to a depth of 15 cm, machine seeding, fertilization, harvesting, and stubble removal.
Experimental plots and fields were seeded with maize, which is the main crop in this region.The applied fertilizer rates were 120 units of N ha −1 for urea and 90 units of P 2 O 5 ha −1 for diamonic phosphate.

Point-wise soil sampling
Ground-based point-wise measurements of four soil properties: the apparent dielectric constant (K a ), volumetric water content (θ i ), gravimetric water content (W i ), and mechanical resistance (γ ) profile, were accomplished following two dif-  ferent sampling strategies.In the experimental fields, sampling was performed each five meters in north-south direction across a rectangular grid (50 m by 15 m), resulting in 44 sampling points for each field.In the seven small plots, all variables were collected each 4 m following the same direction along 22 m length transects, resulting in 5 samples per plot and 35 sampling points for the whole sampling area.The study was accomplished during the dry period in 2001 and repeated in 2004, when the Andosol's moisture was close to the steady-state condition.
Non-destructive sampling techniques were followed for all soil properties except gravimetric water content.For W i measurements, the thermo-gravimetric method was used.This technique requires soil sampling and oven drying (105 • C) (Smith and Mullins, 2001).The relative apparent dielectric constant and volumetric water content were measured by TRASE Time Domain Reflectometry (TDR) equipment (Soil Moisture Equip.Corp.), using a 0.15 m long twoguide probe.TDR is sensitive to the soil's dielectric permittivity and measures the water-filled pore space.Since water has a significantly higher dielectric constant than most solid soil matrices, the overall dielectric constant is highly dependent on the water content and a strong correlation exists between the permittivity values and volumetric water content (Dirksen and Dasberg, 1993).TDR measures the soil's dielectric constant over a broad frequency band (between 100 and 1000 MHz, Topp et al., 1980;Dalton and van Genuchten, 1986).
Soil mechanical resistance was measured, by RIMIC CP 20 ultrasonic cone penetrometer.The measurements were done at 15 mm depth increments from the soil surface to a depth of 600 mm, close to each sampling points where K a and W i were measured.A complete profile consisting of 40 measurements at each sampling point was compressed to a single Hurst exponent value by fractal procedures, characterizing the roughness of the mechanical resistance profile.

GPR measurements
The GPR survey was performed during the dry season, several months after crop harvesting, when the rain stopped and the soil reached a near steady-state moisture profile.In the survey the ZOND-12 equipment (Zond-12c, 2000) was used in continuous reflection mode, with an antenna of 2 GHz central-frequency.The radar system unit is carried across the field manually by operator.The maximum operating depth of this antenna varied in studied agricultural fields from 1.0 to 1.2 m.Once interacting with the subsurface the originally uniform pulse becomes rough, and enriched with multiscale information of distinct nature, including noise (Fig. 2).The GPR traces selected for fractal roughness measurements were extracted from the collected common offset georadargrams at every 2 m in the experimental fields and at every 1 m in the plots.This resulted in 88 and 11 sampling points, respectively.Approximately half of the sampled GPR signals were from locations where K a , W i , and γ had also been sampled.Reference velocities of radar waves and the depth of their penetration were estimated using the permittivity values measured in situ by TDR.

Fractal analysis of the GPR records
We considered the amplitude sequence of the backscattered radar pulse as GPR fractal signature and measured its roughness by the Hurst exponent (H ).The Hurst exponent is related to the graph's fractal dimension by H =2−D (Mandelbrot, 1983).
In order to convert the amplitudes of a radar signal to a fractal elevation profile, we selected the wiggle trace from the georadargram, pre-treated it using Paint Shop Pro (7.04), and rotated it to horizontal position.The obtained image (Fig. 2c) was binarized (algorithm Binar, written by Parrot, 2000) and the algorithm Curve (Parrot, 2004) translated it into a time series of signal amplitudes (Fig. 2d).The roughness of this series was measured by the Rescaled-Range and wavelets techniques of the Benoit software, version 1.3 (SCION Corp., 2003).

Statistical and geostatistical analyses
From a geostatistical view-point the soil properties and fractal parameters extracted from selected GPR traces are irregularly sampled regionalized variables.In their mapping and interpretation the irregular topography of the experimental fields was taken into account by using the elevation map of this mountainous area as base-map for the construction of the variability maps.The mapping protocol is available for all compared experimental areas.The spatial variation of the data was modeled and interpolated by universal kriging (GS + geostatistics software, version 3.1.7,Robertson, 1998;and SURFER, version 6.01).The maps of dielectric permittivity; gravimetric and volumetric water content; mechanical resistance (for each 15 mm layers and for the whole profile), the Rescaled-Range (D R/S ) and wavelet-(D W ) fractal dimensions, were independently constructed.Most of the raw semivarogram were fitted with isotropic exponential semivariogram models, trying to minimize the nugget value ε 2 and maximize R 2 .For some variables, only a linear model gave acceptable fit to the experimental semivariogram.
Two types of D R/S and D W maps were prepared for the experimental fields and plots, one based on 44 and the other on 88 values, respectively.Comparison of these maps with the kriging maps of all other variables (Wi, Qi and γ ) was accomplished first visually, and then by quantifying the distribution of the main clusters across the variability maps by traditional box counting, and by measuring the roughness of maps by Rescaled-Range and wavelet techniques (Benoit 1.3 software, SCION Corp., 2003).Statistical comparison of all studied variables was accomplished by the One-Way ANOVA software (Statgraphics 5, Plus, 2000).

Experimental fields
No statistically significant differences either in gravimetric (W i ), or in volumetric (θ i ), water content were found by traditional methods in the first 15 cm of Andosol in the compared fields (Fig. 3a-h).The mean W i varied from 50.6% under conventional (CT) to 50.3% under zero tillage (ZT).The distribution of water content was more homogeneous in the conventional field.Only one cluster with higher gravimetric water content (69.4%) was delimited between 25 and 35 m along the sampling transect, for both studied fields, where the wetter clusters have been associated with the notably darker color.A similar, but not so well defined humid cluster can be visualized on the volumetric water content map (Fig. 3d), where some new small moist areas appeared.However, the homogeneity of θ i under CT was still notorious in comparison with ZT, where the kriged map shows a major heterogeneity and larger clusters identified as areas of water concentration in more uniform soil moisture conditions.In case of gravimetric water content, only two big clusters with high moisture content are distinguished, while on the volumetric water map three smaller humid areas are visualized (Fig. 3h).The mechanical resistance profiles (Fig. 4a-f) show statistically significant differences between the compared sites at two depths: 1. 10-20 cm and 2. 40-50 cm.
Higher mechanical resistance of Andosol under the CT is related to the compacted microlayer formed as the result of soil disking at the 10-20 cm depth, while the more compacted horizon (at 40-50 cm) coincided with the thin plough pan, typical for this agricultural zone.Both dense microhorizons were observed directly on the soil profile opened beside the experimental field.
On the three variability maps obtained for the discussed soil properties (W i , θ i and γ ) the location of humid clusters always coincides with the areas of higher density.Still, no statistical correlation was found between the different pairs of the above-mentioned soil properties.However our data are in agreement with results documented by Veronese Júnior et al. ( 2006) who have concluded that the compacted zones of the plot show higher water content values.
The Rescaled-Range-, as well as the wavelet fractal dimensions of the GPR signatures (Fig. 5a and b) were significantly different for Andosol under CT and ZT (Fig. 5).Wet and dense clusters always corresponded to higher Hurst exponents (lower roughness) of the GPR traces.Higher D R/S and D W were always more probable for the conventionally managed Andosol.While strikingly different signatures, and very different Hurst exponents, can be observed in waves backscattered from soils with different water content and mechanical resistance (Fig. 6), no significant correlation was found between the measured fractal dimensions and water www.nonlin-processes-geophys.net/15/711/2008/ Nonlin.Processes Geophys., 15, 711-725, 2008  content, or mechanical resistance.We speculate that each fractal dimension encodes at the same time the information about a set of soil properties, so it is not so easy to single out the individual role of any property by the traditional analytical techniques.The higher macrofauna activity in ZT (Oleschko and Chapa, 1989) and therefore the major presence of big cavities inside this agricultural field may also contribute to this lack of data correlation.
Fluctuations in the values of D R/S and D W were significantly different for the compared management practices (Fig. 5).For instance, in CT, the mean value of D R/S was 1.592, decreasing from 1.674 in dry-to 1.513 in wet clusters.In ZT, the same D R/S has a mean of 1.628, ranging from 1.779 to 1.530 in dry, respectively to humid clusters.Note, the values of fractal dimension in dry soil with higher moisture content, extracted from the GPR traces was always closer to the D of white noise.
In order to have a comparable sampling rate between point-wise measured properties and microwaves extracted from the georadargrams, initially we mapped only 44 fractal dimensions for each field.The clusters on the microwavederived D R/S and D W maps coincided visually with those on the water-content and mechanical resistance maps (Compare Fig. 6 with Figs. 3 and 4).It should be mentioned that when Weihermüller et al. (2007) intended to map the spatial variation of soil water content at the field scale with different ground penetrating radar techniques, using the indirect moisture estimation from the wave velocity in porous media, the results were difficult to interpret due to the strong attenuation of the GPR signal.This attenuation was related by these authors to the silt loam texture at the test site.
Notwithstanding, when all fractal dimensions extracted from GPR traces (altogether 88) in the present research, have been utilized for the construction of more detailed variability maps, the position of W i , Q I = i and γ and clusters inside the fields did not change significantly but the areas of different dimensions were better delineated (Fig. 7).
In order to map the measured soil properties and fractal dimensions, with one exception, the raw semivariograms were fitted with isotropic exponential models (Nugget, sill, range and goodness of fit are compiled in Table 1).The semivariogram of the 88 values of the wavelet dimension D W was fitted by a linear model.A monotone increase of the semivariograms with increasing lag (except in case of the D W values, Table 1) and their asymptotic tendency to a constant maximum "sill" have been observed for all compared variables.

Statistical comparison between spatial variability maps
At the next step of this study each map was analyzed as gray scale image (.bmp).The roughness of gray tones distribution across the spatial variability maps of studied Andosol physical and mechanical properties (W i , θ i and γ ) was measured by Hurst exponent extracted by Rescaled-Range and wavelet techniques and compared with the roughness of D R/S and D w maps.Significant statistical correlation was found between all pairs of compared variables except γ 10−20 /D R/S and W i /D w by r-Pearson analysis (Table 3), agreed with no significant differences between all pairs of compared variables shown by t-Student test.Therefore, the location of dark clusters (with high variable values) not only coincided when maps were compared visually but are also characterized by statistically similar gray intensity roughness.

Experimental plots
In the seven small experimental plots the spatial variation of the gravimetric water content W i , dielectric permittivity (K a ), volumetric water content (θ I ), and mechanical  resistance (γ ) at 10-20 cm depth shows a clear correlation with the applied management system (Table 2).However, not all variables present statistically significant differences.The comparisons by One-Way ANOVA (Statgraphics 5, Plus, 2000), proved the advantages of using the fractal dimensions instead of the other variables.When the dielectric permittivity, volumetric water content, and mechanical resistance of Andosol are compared by pairs for seven experimental plots, only the soil under CT was significantly different (with respect) to the rest of the management systems (at 95.0% confidence limit).Unclear, and sometimes counterintuitive, trends were documented for all other plots.On the other hand, the information extracted from the fractal dimensions of the GPR traces was always more precise.It was not only consistent with all trends of soil dynamics documented in Table 2, but also statistically significant for most of the compared management systems (Fig. 8).First of all, both fractal dimensions detected differences in Andosol with and without crop residues roughness.Significantly higher roughness (higher D-and lower H-values) is observed in subsurface structure of plots with zero tillage and 33, 66 and 100% of crop residues.For instance, the GPR traces extracted from ZT-100 plot have 2) * The coefficient of variation (%) is in brackets.ZT-0 = zero tillage without crop residues.ZT-33, ZT-66 y ZT-100 = zero tillage with 33, 66 and 100% crop residues, respectively.MT = minimal tillage CT = conventional tillage.
the highest roughness in comparison with all other management treatments, which have been expressed in lower Hurst exponent value (0.17) in comparison with Andosol under zero tillage and without organic cover (H =0.42).The roughness of waves scattered from the bare (H =0.44) and under conventional tillage (H =0.47) soil was similar.Therefore, both fractal dimensions D R/S and D W carry more useful information about the explored porous media than the other traditional measurements.However, sometimes one of these fractal parameters worked better from statistical view-point than the other.The coefficient of variation was always higher for wavelet fractal dimension, duplicating or sometimes triplicating the value of D R/S coefficient of variation.The value of D R/S (1.56), obtained for the bare plot with higher volumetric water content (26.6%) and larger compaction (1.5 MPa at first 10 cm), was comparable with the D R/S of ZT-0 and CT (1.58 and 1.53, respectively), while the D w values for the same three plots showed significant differences.In this case, the wavelet fractal dimension has the smallest value (and consequently the lowest roughness -largest Hurst exponent) under zero tillage without residues cover, and under conventional management (1.22 and 1.28, respectively).This trend is in agreement over the tendencies of K a -, θ i,and γ -dynamics, being both D w significantly smaller than in the bare plot (where D w is 1.42).Therefore, in some cases the wavelet dimension fits better to the visual field perception of the trends, and it is a clear indicator of the Andosol spatial variability.However, for seven of the 21 compared pairs of plots, differences in the wavelet fractal dimensions were not statistically significant.Curiously, three of these seven pairs were significantly different with respect to the Rescaled-Range fractal dimension.Considering the tendencies observed in nearby plots for the dynamics of both dimensions, there were no significant differences in Andosol physical properties between zero tillage with 33% organic residues vs. zero tillage with 66% organic residues, and zero tillage with 66% vs. zero tillage with 100% cover.Beside the above-mentioned exceptions, the statistical analysis of the GPR trace roughness has detected significant differences in the soils' scattering properties even under very similar tillage conditions.
In the experimental plots, the fractal analysis was applied in order to compress the forty mechanical resistance values of each penetrometer' profile inside the single value of fractal dimension.The range of the D W semivariogram for mechanical resistance (r=21.3m) was similar to the range of the D R/S extracted from the GPR signature (r=22.8m).The ranges for the permittivity (r=15.8m) and gravimetric water content (r=12.9m) semivariograms were much smaller.The largest range (r=37.5 m) corresponds to the semivariogram of volumetric water content, while the lowest (r=9 m) to that of the wavelet fractal dimensions extracted from the GPR traces.The maximum mean value for W i (65.4%) was found in the ZT-33 plot, while the biggest K a (15.8) and θ i (30.0%) corresponds to the ZT-66 field.The highest mechanical resistance γ 50−60 cm =2.5 MPa was documented for the Andosol under zero tillage with 33% of crop residues.This plot, together with the Andosol under zero tillage without crop residues (Z0-0, K a =14.3, θ i =27.9% and γ 0−10 cm =1.3 MPa), have the most homogeneous distribution of mechanical resistance inside the sampled profiles at the depth from 10 to 40 cm.The lowest coefficients of variation were documented for gravimetric (0.8%) and volumetric (6.6%) water content of Andosol under zero tillage without crop residues experimental plot.Notwithstanding, if the variations of all studied parameters are compared, the Rescaled-Range fractal dimen-   sion would be visualized as the less variable: CV is fluctuating between 1.9 and 9.6 for all experimental plots.The variation of wavelet dimension was higher, increasing from 4.7 under zero tillage to 17.7 under zero tillage with 66% of crop residues.The mechanical resistance of Andosol at depth of 0-10, 10-20, 20-30 and 30-40 cm was the parameter with highest coefficients of variation (58.8; 59.5; 62.9 and 65.1%, respectively) under zero tillage with 100% of crop residues.
As mentioned above, all compared variables were sampled regularly along the reference transects across the experimental plots.Therefore, at the beginning we constructed the semivariograms and the kriged maps independently for each plot.However, in Fig. 10, we superimposed all maps in order to enhance the detected tendencies in fractal dimension dynamics under the compared tillage practices.The tendency for homogenization of the studied properties has been documented and discussed above for the field under CT and ZT.Similar trends were observed on the maps constructed either from the volumetric water and mechanical resistance values, or from the fractal dimensions extracted from GPR traces.The highest variability of K a , θ i, and γ was detected under the zero tillage plot with 66% of crop residues (ZT-66, Table 2).The same is seen on the D R/S map extracted from GPR records.
Three plots under zero tillage, with 33, 66, and 100% residues, show strong clustering in the spatial distribution of permittivity, moisture, and mechanical resistance, coinciding with the slope direction of the agricultural fields.Clusters with similar morphology are observed on the D R/S and D W maps too.Under these management conditions the highest volumetric water content was found.The discrete θ i values have fluctuated between 24.6 and 37.1%.High water content (29.2%) occurs in the same area as the maximum mechanical resistance (2.5 MPa), of the lower part of the zero tillage plot with 33% of residues.Visually, this relation appears as a cluster of high-values on the γ -map (for the 0-10 cm layer), and on the map of Rescaled-Range fractal dimension.
For better visualization of the differences between the spatial structures of studied variables, we used for mapping all experimental data coming from the seven compared plots, assuming that they were sampled in a statistically homogeneous region.As, in our case, once all experimental lots used to belong to the same agricultural field under conventional management, the assumption is possibly not far-fetched.At present, no statistically significant differences can be found between the plots in terms of their common physical and mechanical properties.Therefore, we used 77 data to krige and map each studied property.The complete data set fitted well the exponential semivariogram model, with increased values of R 2 .As compared with the above-discussed analyses performed separately for each plot, this study on the complete data set detected spatial variations of the compared properties within shorter sampling intervals (Fig. 10).
The smallest values of Rescaled-Range fractal dimensions, and therefore, the highest Hurst exponents and smallest microwave roughness, were found for the bare Andosol, for the soil under zero tillage without residues, and for the soil under conventional tillage (1.56, 1.58 and 1.53).The highest roughness was typical for microwaves scattered from the Andosol under zero tillage with 100% organic cover.
Note that the tendencies in Andosol dynamics under tillage documented for small experimental plots were similar to those documented for the two big experimental fields and discussed above.

Conclusions
We documented that the fractal dimension of EM microwaves backscattered from soil is dependent on, and strongly correlates with, the heterogeneity of the soil's physical properties.Our proposed approach for non-invasive soil mapping provides a fast, cheap, and non-destructive way to delimit the coherent clusters and structural patterns of the soil's spatial variability.In this research we checked the FDM technique on different plots of volcanic soil (Andosol) under contrasting management practices.The differences in measured permittivity, volumetric and gravimetric water contents and mechanical resistance were not always statistically significant between soils under different management systems.In seven experimental plots, the comparison of these variables by the ANOVA statistical analysis detected differences only between the conventional tillage as compared with each of the other six management systems.However, when the fractal dimensions of the recorded GPR traces were mapped, the slight changes in soil properties were en-hanced and the differences have become statistically significant.The fractal dimension maps show significant differences between all pairs of the compared tillage systems, except three of them most close in design to each other.These results show the advantages of the proposed, GPR-based, continuous Fractal Dimension Mapping of soil variability, and the ability of this technique to detect small changes in the soil's physical and mechanical properties.We hope that the technique would open up new vistas for the multiscale non-destructive mapping of soil variability, and for more extensive GPR use by Soil Service Agencies (Doolittle, 1987;Dwivedi, 2001).

Fig. 3 .
Fig. 3. Semivariograms of the Andosol's gravimetric (a and e) and volumetric (c and g) water content in the Conventional Tillage (a and c) and Zero Tillage (e and g) fields.Clusters with different water content appear on the maps of soil variation under CT (b and d) and ZT (f and h).All maps, constructed by universal kriging, are superimposed on the topographic map.

Fig. 6 .
Fig. 6.Spatial variability maps (a and d) and the best-fitting exponential semivariogram models (b and e) of the Andosol's GPR-derived Rescaled Range (D R/S ) and wavelet (D W ) fractal dimensions, constructed from 44 sampling points for CT and ZT fields, respectively.Some examples (c and f) of GPR traces used for fractal analysis.The D R/S (a) and D W (c) clustering is comparable with the soil property maps (Figs. 3 and 4).

Fig. 7 .
Fig. 7. Best-fitting exponential semivariogram models (a and c), spatial variability maps (b and e) of the Andosol's rescaled range (D R/S ) and wavelet (D W ) dimensions in the CT and ZT fields.Altogether 88 traces were extracted, 2 m apart.The D R/S and D W (a anf c) clustering can be compared with the maps of soil physical properties (Figs. 3 and 4).

Fig. 9 .
Fig. 9. Spatial variability map of the Rescaled Range fractal dimension D R/S along the seven compared plots.For each plot the map was constructed by universal kriging, independently of the others.All maps were combined to make clear the fractal dimensions' clustering in the differently managed plots.(Management codes are explained in Table 2.)

Fig. 10 .
Fig. 10.Exponential models for the Rescaled-Range (a) and wavelet (d) fractal dimension semivariograms, constructed from all data.The spatial variability maps for the GPR traces' D R/S (b) and D W (e) dimensions were calculated by universal kriging, using all experimental data together.Observe the clear difference between GPR traces recorded over different plots(c).

Table 1 .
Exponential model parameters for the semivariograms over the experimental fields.

Table 2 .
Mean values of the variables in Andosol under seven management systems.

Table 3 .
r-Pearson analysis of spatial variability maps roughness for gravimetric (W i ) and volumetric (Q i ) water content, mechanical resistance (γ (10−20 cm) ) vs. roughness of D R/S and D w maps.