Multifractal analysis of 3D images of tillage soil

ARTICLE IN PRESS GEODER-12630; No of Pages 8 Geoderma xxx (2017) xxx–xxx

Contents lists available at ScienceDirect

Geoderma journal homepage: www.elsevier.com/locate/geoderma

Multifractal analysis of 3D images of tillage soil I.G. Torre a, b, * , J.C. Losada b , R.J. Heck c , A.M. Tarquis b, d, e a

Department of Applied Mathematics and Statistics, EIAE, Polytechnic University of Madrid (UPM), Madrid, Spain Complex Systems Group, UPM, Ciudad Universitaria sn, Madrid 28040, Spain School of Environmental Sciences, Ontario Agricultural College, University of Guelph, Canada d Department of Applied Mathematics, ETS Agronomic Engineers, Polytechnic University of Madrid (UPM), Avenida Puerta de Hierro 2, Madrid E28040, Spain e Research Centre for the Management of Agricultural and Environmental Risks (CEIGRAM), Polytechnic University of Madrid (UPM), Calle Senda del Rey 13, Madrid E28040, Spain b c

A R T I C L E

I N F O

Article history: Received 7 April 2016 Received in revised form 7 January 2017 Accepted 16 February 2017 Available online xxxx Keywords: Multifractal 3D-images Gliding method Cube Counting Spatial distribution

A B S T R A C T Recently, X-ray microtomography (lCT) has open a new way to study soil pore structures. However, lCT data originally comes as gray scale images and the selection of segmentation method to binarized it has an important influence in the structure characterization. Three soil lCT 3D images, corresponding to ploughed soil with different tillage tools, were used in this study. A multifractal analysis was applied to the original gray images avoiding image binarization to characterize and differentiate different soil structures. In this analysis we took into account the effect of image resolution and the subdivision method applied, Cube Counting (CC) and Gliding Cube (GC). Comparisons among the multifractal spectrums (MFS) estimated indicated that the reduction in resolution affected more in Molboard image in which the minimum scale included in the analysis was one voxel and it was almost imperceptible in the Chisel image. The three multifractal spectrums were quite distinctive when the maximum resolution was chosen and the GC method was used, meanwhile the CC method obtained the MFS of Chisel and Roller too close to be differentiated. The amplitude and symmetry of the multifractal spectrums point out the influence of each tillage tool in the hierarchy of soil structure. Moldboard creates a higher complexity in soil structure as it physically removes the soil. Chisel tends to destruct the aggregates homogenizing soil structure and presents a weak multifractal nature. Finally, Roller is an intermediate case with a scaling character mainly in the lower gray values of the soil image. © 2017 Elsevier B.V. All rights reserved.

1. Introduction There is increasing evidence that quantitative characterization of the soil structure and of its heterogeneity holds the key to a deeper understanding on physical, chemical and biological processes that take place within them (Young and Crawford, 2004; Blair et al., 2007). Recent advances in imaging techniques, such as X-ray Computer Tomography (X-ray CT) have opened new perspective for the quantification of the internal pore structure of soils as well as characterizing their heterogeneity (Pierret et al., 2002; Anderson et al., 2003; Rachman et al., 2005; Gibson et al., 2006b; Taina et al., 2008; Garbout et al., 2012). X-ray CT provides qualitative and quantitative information on the 3D physical structure and the spatial distribution of the pore space (Houston et al., 2013) and reveals the extraordinary complexity of this system (Santiago et al., 2008; Cárdenas et al., 2010). At the same time, these techniques * Corresponding author at: Department of Applied Mathematics and Statistics, EIAE, Polytechnic University of Madrid (UPM), Madrid, Spain. E-mail address: [email protected] (I. Torre).

reduce the physical impact of sampling and allow to a rapid scanning to study sample dynamics in near real-time (Garbout et al., 2013a). This is a significant step to study natural porous media at micro-scale. Therefore, X-ray CT soil images can give a good understanding of soil structural changes due to soil tillage and compaction (Taina et al., 2008). We can found in the literature several studies documenting these soil management effects on soil physical properties using X-ray CT systems (Olsen and Børresen, 1997; Munkholm et al., 2003; Papadopoulos et al., 2009; Garbout et al., 2013b; Munkholm et al., 2013 among others). For example, Gantzer and Anderson (2002) highlighted the number and size of macropores importance in tillage-induced structures. Rasiah and Alymore (1998) found that macropore physical properties affect water flow rate and retention using an X-ray CT analysis. Fractal geometry has been increasingly applied to quantify soil structure, using fractal parameters derived directly through image analysis, due to the complexity of the soil structure, and thanks to the advances in computer technology (Tarquis et al., 2003). First, several works have been done extracting mass fractal and surface

http://dx.doi.org/10.1016/j.geoderma.2017.02.013 0016-7061/© 2017 Elsevier B.V. All rights reserved.

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013

ARTICLE IN PRESS 2

I. Torre et al. / Geoderma xxx (2017) xxx–xxx

fractal dimensions (Brakensiek et al., 1992; Pachepsky et al., 1996; Giménez et al., 1997; Gimenez et al., 1998; Oleschko et al., 1998; Bartoli et al., 1999; Dathe et al., 2001; Dathe and Thullner, 2005; Rogasik et al., 1999; Gantzer and Anderson, 2002; Perret et al., 2003). There are many studies on soil tillage research applying fractal methods focused on macroporosity or pore surface area based on CT-scan soil images (Perfect and Kay, 1995 see references therein). Later on, the studies based on soil images have evolved to a multifractal analysis -MFA (Posadas et al., 2003; Bird et al., 2006; Dathe et al., 2006; Martínez et al., 2010) becoming widely applied. Even we focused this work on the use of MFA to study the relation of the characterization of soil structure among scales, there are several works studying wider scaling behaviours. These analyses include extended power-law scaling (linear relations between log structure functions of successive orders) at all lags, and frequency distributions of the variable’s increments tending to be symmetric with peaks that grow sharper, and tails that become heavier, as the lags between pairs of values decreases (Guadagnini et al., 2014). However, searching on multifractal methods in soil tillage research almost of them are dedicated to soil microtopography (Vázquez et al., 2008), soil surface roughness (Roisin, 2007; García Moreno et al., 2008) and quantify distributions of soil properties (Siqueira et al., 2013; Jonard et al., 2013). Despite these progresses, there remains a lack of general agreement to the appropriate pore-solid CT threshold (Tarquis et al., 2008; Cortina-Januchs et al., 2011), which is used to obtain a black and white image from the original gray scale data, before calculating any of these parameters. Gibson et al. (2006a) compared three fractal analytical methods to quantify the heterogeneity within soil aggregates; in this work, the frequency distribution of pore and solid components was clearly dependent on thresholding, which could not be generalized. Tarquis et al. (2009) point out that a practical problem in the MFA of binary images is that the thresholding method have a pronounced effect on the porosity and resulting generalized dimensions. It has been suggested to further study grayscale soil images for multifractal characterization of soil structure avoiding any intermediate thresholding step (Zhou et al., 2010; Zhou et al., 2011). This MFA applied in grayscale images is very common in other type of images (Lovejoy et al., 2008; Tarquis et al., 2014). On the other hand, MFA normally involves partitioning the space of study into non overlapping boxes (cubes in the case of this study) to construct samples with multiple scales, known as the box counting method and named in this study as Cube Counting (CC) method. The number of samples at a given scale, applying it, is restricted by the size of the partitioning space and data resolution, which is usually another main factor influencing statistical estimation in MFA (Cheng and Agterberg, 1996). To avoid these problems several methods have been proposed being one of them the Gliding box method (GB) already applied in several type of images (Grau et al., 2006; Tarquis et al., 2007) and in this study we will named it as Cube Gliding (CG) method.. The two objectives of this work are to study: 1) the effect of the subdivision method applied in the calculation of the multifractal parameters, comparing the CC and CG method and the effect of image resolution; 2) the effect of the different soil structure in the multifractal parameters calculated on gray scale soil images. With these purposes, 3D CT scan gray scale soil images, extracted from areas ploughed with different tools, were used.

The trials were conducted at the Alameda del Obispo experimental farm (38N, 5W, altitude 110 m), Cordoba, Spain (see Fig. 1). The climate is Mediterranean with a mean annual rainfall of 595 mm. Summer in Cordoba is dry and hot while autumn and winter are mild and rainy (Peel et al., 2007). The soil is a loamy alluvial with particleg size distribution in the upper (0–15 cm) soil layer: sand, 350 kg ; silt, g g 443 kg ; and clay 207 kg (Boulal et al., 2011). It is classified as Typic Xerouvent (Soil Survey Staff, 2010). Intact cores of soil samples for each of the tillage treatment were packed into polypropylene cylinders of 8 cm diameter and 10 cm high. These were imaged using an mSIMCT at 155 keV and 25 mA. An aluminium filter (0.25 mm) was applied to reduce beam hardening and later several corrections where applied during reconstruction. All 3D volumes were converted using VGStudioMax v.1.2.1 into image stacks with voxel-thick slices. All soil samples were scanned and reconstructed into 3-D volumes with a voxel size of 80 lm. From the reconstructed volumes a centred volume of 256 × 256 × 256 voxels (20.5mm × 20.5mm × 20.5 mm) were selected for all samples. Fig. 2 shows the original 3-Dimensional model reconstructed from the slices obtained with CT-scan. We have selected a 256 × 256 × 256 voxel subsample for the analysis of each different tillage. The full size is computed through the Cube Counting and Gliding method, and then scaling it with Bilinear Interpolation algorithm using ImageJ 3D software to 128 × 128 × 128 voxel and analysing again with both methods. 2.2. Multifractal analysis An object has its topological measure which is one of the several ways of defining the dimension of the space. For example, the topological dimension of a line is one,the dimension of surface is two, and the dimension of a 3D object is 3. However there are some natural and mathematical objects that exceeds their topological dimension. A geometrical fractal is a mathematical object which can exceed its topological dimension and typically displays self-similar patterns, e.g., the borders of a country (Mandelbrot, 1967), the Koch snowflake (Mandelbrot, 1983) or Sierpinski triangle. When we pass from a geometrical object to a measure, as gray values in an image could be interpreted, we are passing to another type of analysis that can be seen as an extension of fractals. Multifractal analysis initially appeared to study of energy dissipation on multiplicative cascades models in the context of the fully developed turbulence (Mandelbrot, 1999). After that it has been implemented on several different natural systems.

2. Material and methods 2.1. Study area and image acquisition For this study we have analysed 3 real soil samples. All of them have been extracted from the same experimental farm in Cordoba (Spain) but ploughed with different tillage tools.

Fig. 1. Soil samples were extracted from different areas ploughed with three different tools: Moldboard in yellow, Chisel in green and Roller in blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013

ARTICLE IN PRESS I. Torre et al. / Geoderma xxx (2017) xxx–xxx

(a) Moldboard tillage

3

(b) Roller tillage

(c) Chisel tillage

Fig. 2. Visualization of the soil samples used in this study. Each volume is build with 256 × 256 × 256voxels and each voxel size is 80lm.

There are several methods to subdivide the set in which we are going to compute the fractal dimension of an image; each one has its own theoretical basis. Here, we are going to use two of them: Cube Counting (CC) and Cube Gliding (CG), which is an improved method for the low density areas of the image (Saa et al., 2007). The Cube Counting methodology defined by Russell et al. (1980) for box-counting, is a classical in this field. It breaks the spatially extended dataset into smaller and smaller pieces (cubes) and analyse them at each smaller scales. One way to understand this method is as zoom in to observe how details change with scale. In the algorithm implemented, the cubes are scaled in power of two. Therefore, we use images which its size in voxels is a power of 2, one will be 256 and other one 128. In the multifractal method, a weighted factor (mass exponent, q) is applied to datasets extracted from the image giving more or less importance to the high or low mass density areas of the picture. Meshing the 3D image with cubes of size 4 and define Pi (4) as the probability of this mass at the cube i relative to the total mass of the picture for a cube size 4, then the generalized dimensions Dq which corresponds to scaling exponents for the qth moments of the measure are defined as (Meneveau and Sreenivasan, 1987):

Dq =

log

1 lim q − 1 4→0



q

i Pi (4)

and defining the partition function as X(q, 4) = ized dimensions can be expressed as:

Dq =

(1)

log 4

log X(q, 4) log 4



q i Pi (4)

the general-

(2)

for q = 1. For q = 1, D1 is defined by taking the limit when q approaches 1 and by using L’Hoplital’s Rule (Teng et al., 2010):  D1 = lim

i Pi (4) log Pi (4)

4→0

log 4

(3)

f (a) = qa − t

(5)

dt dq

(6)

a=

An easily used definition of the singularity spectrum which is the one used is (Chhabra and Jensen, 1989):  f (q) = lim

i li (q, 4) log li (q, 4)

 a(q) = lim

(4)

for q = 1 and in the case of conservative measure t1 = 0. Another interesting measure in multifractal analysis is the relationship between a Hausdorrf dimension f and an average singularity

(7)

log 4

4→0

i li (q, 4) log Pi (q, 4)

(8)

log 4

4→0

where l i (q, 4) is how the distorted mass probability at a cube compares to the distorted sum over all cubes at this cube size: (Pj (4))q li (q, 4) =  q j (Pj (4))

(9)

The disadvantage of the CC method is that for negative q values the uncertainties are high and the errors in the estimation are bigger. CG method has been developed to avoid this problem. It was originally used in lacunarity analysis and later modified by Cheng in (Cheng, 1999) to apply it in the multifractal analysis. Later, it was used in soil binarized image analysis (Grau et al., 2006). This method constructs samples by gliding a cube of certain dimensionless size 4 over the grid map in all possible directions creating overlapping cubes. CG estimates first the mass exponent function t(q): log t(q) + D = lim

4→0

This measure Dq is related with the mass exponent t(q) as follows: tq = (q − 1)Dq

strength a as implicit functions of the parameter q. Those are defined as follows:



1 N∗ (4)

N∗ (4) i=1

q



li (4)

log 4

(10)

where 4 is the dimensionless cube size, N∗ (4) represents total number of gliding cubes of size 4 with measure l(4) = 0 and D is the topological dimension of the object analysed (D = 2 for two-dimensional objects and D = 3 for three-dimensional objects) The functions a and f(a) are estimated numerically by the Legendre transformation as a(q) =

dt(q) , dq

f (a) = a(q)q − t(q)

(11)

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013

ARTICLE IN PRESS 4

I. Torre et al. / Geoderma xxx (2017) xxx–xxx

We have calculated the error of a knowing the uncertainty of t and the explicit relationship between the two variables (Taylor, 1997). The advantage of gliding cube method is that the larger sample size provides better statistical results. The disadvantage is that it is more computationally demanding due to the higher number of cubes that generates and where calculations on the measure has to be applied. In this study we have applied both subdivision methods (CC and GC). The range of q values used has been from q = −10 to q = 10. The scales chosen for the linear fit were selected observing the linearity of the bilog plot of the partition function so error was minimized. 3. Results and discussion 3.1. Descriptive statistics on the gray values The histograms of the gray values for each tillage tool is showed in Fig. 3 in which Moldboard presents a clear bimodal shape meanwhile Roller and Chisel histograms are more symmetrical and they do not present two local maximum. In the case that we wanted to performance a binarization of the images applying a global thresholding, based on the histogram, it will be very difficult in the former cases. Respect to the moments of the distribution of gray values (Table 1), again Moldboard presents a higher standard deviation and there is a significant difference between mean and mode indicating its asymmetry. On the other hand, Roller and Chisel have lower standard deviation and their mean and mode are practically the same. On the other hand, Roller and Chisel have lower variance and their mean and mode are practically the same corresponding to almost null skewness. Comparing the kurtosis values from each tillage tool, Roller has the highest value indicating a tail shape significantly different from a normal or Gaussian distribution. 3.2. Comparison between subdivision methods and resolutions The first important step in the multifractal analysis is to determinate which range of scales is going to be used in the calculations (Saa et al., 2007). This depends mainly in the linearity that the partition function (w(q, 4)) presents in the bi-log plot (see Fig. 4). For all of the cases the largest scales are selected, however for each tool different minimum scale have to be chosen. In Molboard case the minimum scale is 1 voxel, Roller presents a clear linearity till 8 voxels and Chisel for 4 voxels of cube length size. The dashed lines in Fig. 4 mark the smallest scale considering to fit a linear regression. These selected ranges were used for both methods, Cube Counting (BC) and Gliding cube (CG). In the CC method the multifractal spectrum (MFS) was calculated using Eqs. (7) and 8 and the standard errors of the estimations are

Fig. 4. Analysis of w(q, 4) vs log(4) to determine the appropriate range of scales selected for the multifractal analysis. From top to bottom the minimum size selected (marked by dashed line) has been: 1 voxel for Moldboard, 8 voxels for Roller and 4 voxels for Chisel.

coming from the errors in the slopes. In the GC method the spectrum is derived from t(q) using Eq. (11) and the errors are coming from the uncertainty of the mass function. In the case of Moldboard, the right part of the MFS, that corresponds to negative q values, change between the two methodshas been previously reported (Grau et al., 2006). This difference is more evident at 128 × 128 × 128 voxels resolution than at 256 × 256 × 256. This is mainly due to the change in linearity of the bilog plot of w(q, 4), when the scales are selected for the linear regression. However, the shape of the multifractal spectrum is similar in all cases showing an asymmetry in its left part. This is pointing to a stronger scaling behaviour in the high gray values than in the low ones; these last ones correspond to the darker voxels. Respect to Roller, CC method estimates a MFS almost identical at both resolutions and quite symmetrical. The range of values achieved

Fig. 3. Histogram of the gray value of each 3D image: Moldboard, Roller and Chisel. Higher values of pixel intensity indicate lighter gray and lower values are for darker pixels.

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013

ARTICLE IN PRESS I. Torre et al. / Geoderma xxx (2017) xxx–xxx Table 1 Statistical moments of the gray values distribution of each 3D image: mean, mode, variance, kurtosis and asymmetry. The confidence intervals for the mean and variance are stated at the 95% confidence level. Levene test for the variance and Welch’s t-test for the means have been passed showing statistical differences across samples in both cases with p-values< 0.01. Sample

Moldboard

Roller

Chisel

Mode Mean Variance Skewness Kurtosis

119 112.251 ± 0.006 151.6 ± 0.1 −0.930 0.309

102 101.219 ± 0.003 70.07 ± 0.04 0.003 2.415

101 101.096 ± 0.003 58.44 ± 0.03 −0.021 0.964

in a and f(a) is much more reduced compared to Moldboard. In the case of CG method there is an important difference in the right part of the multifractal spectrum when the estimations are performance at the two resolutions. With the maximum resolution (256 × 256 × 256 voxels) there is a stronger scaling behaviour at the negative q values showing an asymmetry in the right part of the MFS. At lower resolution (128 × 128 × 128) the MFS is quite symmetric showing both parts similar behaviour. In this case there is a clear influence of the Bilinear Interpolation implemented in ImageJ.

5

Finally, Chisel presents a similar case than Roller (Fig. 5). The MFS obtained by both subdivision methods and at both resolutions are very similar. The multifractal character of the spatial distribution of these gray values is very weak as the spectrum is quite reduced compared with Moldboard and Roller. 3.3. Comparison between the tillage tools Using the images with maximum resolution, we are going to study the differences that can be found through each subdivision method. In the case of CC (Fig. 6) it is quite obvious that Moldboard presents a higher multifractal character that the other two expressing a higher hierarchy in the spatial distribution of the gray values. Roller and Chisel are quite similar, appreciating a small difference between them. Several parameters from the MFS are obtained (Table 2). The amplitude of the multifractal spectrums (a max − a min ) confirms the first visual comparison between Roller and Chisel, there isn’t a significant difference. The symmetry of the spectrum can be quantified through (f(a max ) − f(a min )) in which Roller and Chisel present a negative difference with a value close to null pointing out the same scaling behaviour for higher and lower gray values. Again, Chisel shows an opposite pattern with a positive difference due to

Fig. 5. Multifractal spectrum estimated with two different methods: Cube Counting (blue box), Cube gliding (red circles) at two different resolutions (256 × 256 × 256 in the first column and 128 × 128 × 128 in the second one) for each tillage tool (Moldboard, Roller and Chisel in the first, second and third rows.) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013

ARTICLE IN PRESS 6

I. Torre et al. / Geoderma xxx (2017) xxx–xxx

Fig. 6. Multifractal spectrum of Moldboard, Roller and Chisel tillage tools using Cube Counting subdivision method at 256 × 256 × 256voxels resolution.

Fig. 7. Multifractal spectrum of Moldboard, Roller and Chisel tillage tools using Cube Gliding subdivision method at 256 × 256 × 256voxels resolution.

present a clear asymmetry in its MFS and, therefore, a strong scaling in the high gray values. As we can see, the generalized dimensions calculated for q = 1 and q = 2 (D(q = 1), D(q = 2)) are very similar and in this case cannot help to differentiate among the structures studied here. It is mainly in the extreme values where the differences could be significant and for this reason is important to obtain these values with the minimum possible error. The D(q = 0) value is the Euclidean dimension for the three of them as the support of the measure is not fractal. In the case of CG method (Fig. 7), similar results are obtained. Moldboard shows a wider spectrum that markedly differentiated it from the other two, even the errors in f(a) are higher than in the BC method. This is due to the proportion of the error coming from the

mass function. Roller and Chisel present now a different multifractal spectrum mainly in the right side showing the former a stronger scaling in the low gray values. The comparison of the three MF spectrum shows a distinction among the three tillage tools being the a(q = 0) with a different value in each one of them. Quantifying the amplitude and symmetry of the MF spectrum (Table 3) confirms the graphic result showed earlier. In this case the amplitude of the spectrums (a max −a min ) is significantly different among the three of them. The symmetry of the spectrum (f(a max ) − f(a min )) keeps the sign that showed in BC method making again a difference between Roller and Chisel. Respect to the generalized dimensions calculated for q = 1 and q = 2 (D(q = 1), D(q = 2)) using the CG subdivision method, the results are very similar than the ones using CC method.

Table 2 Extreme values and parameters derived from the multifractal spectrums using Cube Counting subdivision method at a resolution of 256 × 256 × 256voxels for each tillage tool: M (Moldboard), R (Roller), C (Chisel). Support dimension (D(q = 0)), entropy dimension (D(q = 1)) and correlation dimension (D(q = 2)) are showed.

Table 3 Extreme values and parameters derived from the multifractal spectrums using Cube Gliding subdivision method at a resolution of 256 × 256 × 256voxels for each tillage tool: M (Moldboard), R (Roller), C (Chisel). Support dimension (D(q = 0)), entropy dimension (D(q = 1)) and correlation dimension (D(q = 2)) are showed.

Method and resolution

Cube Counting 256 × 256 × 256

Method and resolution

Gliding 256 × 256 × 256

Tillage tool

M

R

C

Tillage tool

M

R

C

a min (q = 10) f(a min ) a max (q = −10) f(a max ) a max − a min f(a max ) − f(a min ) D(q = 0) D(q = 1) D(q = 2)

2.983 ± 0.001 2.909 ± 0.007 3.010 ± 0.001 2.959 ± 0.003 0.026 ± 0.002 0.050 ± 0.010 3.0 2.99930 2.99857

2.997 ± 0.001 2.984 ± 0.001 3.004 ± 0.001 2.979 ± 0.004 0.007 ± 0.001 −0.005 ± 0.005 3.0 2.99984 2.99967

2.997 ± 0.001 2.986 ± 0.001 3.003 ± 0.001 2.984 ± 0.002 0.006 ± 0.001 −0.001 ± 0.003 3.0 2.99986 2.99973

a min (q = 10) f(a min ) a max (q = −10) f(a max ) a max − a min f(a max ) − f(a min ) D(q = 0) D(q = 1) D(q = 2)

2.986 ± 0.001 2.909 ± 0.012 3.014 ± 0.001 2.951 ± 0.016 0.028 ± 0.001 0.042 ± 0.028 3.0 2.99951 2.99897

2.997 ± 0.001 2.982 ± 0.006 3.010 ± 0.001 2.937 ± 0.006 0.012 ± 0.001 −0.045 ± 0.012 3.0 2.99987 2.99973

2.998 ± 0.001 2.987 ± 0.001 3.003 ± 0.001 2.983 ± 0.002 0.006 ± 0.001 −0.004 ± 0.002 3.0 2.99989 2.99978

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013

ARTICLE IN PRESS I. Torre et al. / Geoderma xxx (2017) xxx–xxx

4. Conclusions Advances in imaging techniques, such as X-ray Computer Tomography (X-ray CT) give an opportunity to quantify internal soil pore structure and characterize their heterogeneity and complexity. In this work we have selected three soil images coming from different tillage treatments: Moldboard, Roller and Chisel. A MFA have been applied to study if it is possible to differentiate among the three distinctive structures that each tool creates in the soil based on the original images and avoiding binarization. We have considered the influence of the image resolution and the subdivision method in this type of analysis looking for the best procedure that can better distinct the soil structures. The comparison of multifractal spectrum using cube-counting method shows that Chisel and Roller cases scale similarly. Moldboard presents a different spectrum and a stronger MF character. On the other hand, using Gliding cube method shows almost identical results at 128 × 128 × 128 resolution and at a higher resolution it differentiates better the three tillage tools. Respect to the influence of each tillage tool in the hierarchy of soil structure, Moldboard creates a higher complexity in the soil aggregates as physically removes the soil. Chisel tends to destruct the aggregates homogenizing soil structure and presents a weak multifractal nature. Finally, Roller is an intermediate case with a scaling character mainly in the lower gray values of the soil image. Further research will be conducted to repeat this analysis with several samples and study the consistency of these differences found in the MFS among the three tillage tools. Acknowledgments First author acknowledge the fellowship received by Fundacion ONCE. This research was funded in part by Spanish Ministerio de Ciencia e Innovación (MICINN) through project no. PCIN-2014080, project no. MTM2012-39101-C02-01, and project no.MTM201563914-P.

References Anderson, S., Wang, H., Peyton, R., Gantzer, C., 2003. Estimation of porosity and hydraulic conductivity from X-ray CT-measured solute breakthrough. Geol. Soc. Lond. Spec. Publ. 215, 135–149. Bartoli, F., Bird, N., Gomendy, V., Vivier, H., Niquet, S., 1999. The relation between silty soil structures and their mercury porosimetry curve counterparts: fractals and percolation. Eur. J. Soil Sci. 50, 9–22. Bird, N., Díaz, M.C., Saa, A., Tarquis, A.M., 2006. Fractal and multifractal analysis of pore-scale images of soil. J. Hydrol. 322, 211–219. Blair, J., Falconer, R.E., Milne, A., Young, I.M., Crawford, J.W., 2007. Modeling three-dimensional microstructure in heterogeneous media. Soil Sci. Soc. Am. J. 71, 1807–1812. Boulal, H., Gómez-Macpherson, H., Gómez, J.A., Mateos, L., 2011. Effect of soil management and traffic on soil erosion in irrigated annual crops. Soil Tillage Res. 115, 62–70. Brakensiek, D., Rawls, W., Logsdon, S., Edwards, W., 1992. Fractal description of macroporosity. Soil Sci. Soc. Am. J. 56, 1721–1723. Cárdenas, J., Santiago, A., Tarquis, A., Losada, J., Borondo, F., Benito, R., 2010. Soil porous system as heterogeneous complex network. Geoderma 160, 13–21. Cheng, Q., 1999. The gliding box method for multifractal modeling. Comput. Geosci. 25, 1073–1079. Cheng, Q., Agterberg, F., 1996. Comparison between two types of multifractal modeling. Math. Geol. 28, 1001–1015. Chhabra, A., Jensen, R.V., 1989. Direct determination of the f (a) singularity spectrum. Phys. Rev. Lett. 62, 1327–1330. Cortina-Januchs, M., Quintanilla-Dominguez, J., Vega-Corona, A., Tarquis, A., Andina, D., 2011. Detection of pore space in CT soil images using artificial neural networks. Biogeosciences 8, 279–288. Dathe, A., Eins, S., Niemeyer, J., Gerold, G., 2001. The surface fractal dimension of the soil-pore interface as measured by image analysis. Geoderma 103, 203–229. Dathe, A., Tarquis, A.M., Perrier, E., 2006. Multifractal analysis of the pore-and solid-phases in binary two-dimensional images of natural porous structures. Geoderma 134, 318–326. Dathe, A., Thullner, M., 2005. The relationship between fractal properties of solid matrix and pore space in porous media. Geoderma 129, 279–290.

7

Gantzer, C.J., Anderson, S.H., 2002. Computed tomographic measurement of macroporosity in chisel-disk and no-tillage seedbeds. Soil Tillage Res. 64, 101–111. Garbout, A., Munkholm, L.J., Hansen, S.B., 2013a. Temporal dynamics for soil aggregates determined using X-ray CT scanning. Geoderma 204, 15–22. Garbout, A., Munkholm, L.J., Hansen, S.B., 2013b. Tillage effects on topsoil structural quality assessed using X-ray CT, soil cores and visual soil evaluation. Soil Tillage Res. 128, 104–109. Garbout, A., Munkholm, L.J., Hansen, S.B., Petersen, B.M., Munk, O.L., Pajor, R., 2012. The use of PET/CT scanning technique for 3d visualization and quantification of real-time soil/plant interactions. Plant and soil 352, 113–127. García Moreno, R., Díaz Álvarez, M., Saa Requejo, A., Tarquis, A., 2008. Soil surface roughness analyzed as a multifractal measure. Vadose Zone J. 7, 512–520. Gibson, J.R., Lin, H., Bruns, M.A., 2006a. A comparison of fractal analytical methods on 2-and 3-dimensional computed tomographic scans of soil aggregates. Geoderma 134, 335–348. Gibson, S., Fan, Y., Török, T., Kliem, B., 2006b. The evolving sigmoid: evidence for magnetic flux ropes in the corona before, during, and after CMES. Space Sci. Rev. 124, 131–144. Gimenez, D., Allmaras, R., Huggins, D., Nater, E., 1998. Mass, surface, and fragmentation fractal dimensions of soil fragments produced by tillage. Geoderma 86, 261–278. Giménez, D., Allmaras, R., Nater, E., Huggins, D., 1997. Fractal dimensions for volume and surface of interaggregate pores scale effects. Geoderma 77, 19–38. Grau, J., Méndez, V., Tarquis, A., Díaz, M., Saa, A., 2006. Comparison of gliding box and box-counting methods in soil image analysis. Geoderma 134, 349–359. Guadagnini, A., Blunt, M., Riva, M., Bijeljic, B., 2014. Statistical scaling of geometric characteristics in millimeter scale natural porous media. Transp. Porous Media 101, 465–475. Houston, A., Schmidt, S., Tarquis, A., Otten, W., Baveye, P.C., Hapca, S.M., 2013. Effect of scanning and image reconstruction settings in X-ray computed microtomography on quality and segmentation of 3d soil images. Geoderma 207, 154–165. Jonard, F., Mahmoudzadeh, M., Roisin, C., Weihermüller, L., André, F., Minet, J., Vereecken, H., Lambot, S., 2013. Characterization of tillage effects on the spatial variation of soil properties using ground-penetrating radar and electromagnetic induction. Geoderma 207, 310–322. Lovejoy, S., Tarquis, A., Gaonac’h, H., Schertzer, D., 2008. Single-and multiscale remote sensing techniques, multifractals, and modis-derived vegetation and soil moisture. Vadose Zone J. 7, 533–546. Mandelbrot, B.B., 1967. How long is the coast of britain. Science 156, 636–638. Mandelbrot, B.B., 1983. The Fractal Geometry of Nature. Macmillan Ed.. Mandelbrot, B.B., 1999. Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier. . Multifractals and 1/ NoiseSpringer., pp. 317–357. Martínez, F.S.J., Martín, M., Caniego, F., Tuller, M., Guber, A., Pachepsky, Y., GarcíaGutiérrez, C., 2010. Multifractal analysis of discretized X-ray CT images for the characterization of soil macropore structures. Geoderma 156, 32–42. Meneveau, C., Sreenivasan, K., 1987. Simple multifractal cascade model for fully developed turbulence. Phys. Rev. Lett. 59, 1424. Munkholm, L.J., Heck, R.J., Deen, B., 2013. Long-term rotation and tillage effects on soil structure and crop yield. Soil Tillage Res. 127, 85–91. Munkholm, L.J., Schjønning, P., Rasmussen, K.J., Tanderup, K., 2003. Spatial and temporal effects of direct drilling on soil structure in the seedling environment. Soil Tillage Res. 71, 163–173. Oleschko, K., Brambila, F., Aceff, F., Mora, L.P., 1998. From fractal analysis along a line to fractals on the plane. Soil Tillage Res. 45, 389–406. Olsen, P.A., Børresen, T., 1997. Measuring differences in soil properties in soils with different cultivation practices using computer tomography. Soil Tillage Res. 44, 1–12. Pachepsky, Y., Yakovchenko, V., Rabenhorst, M., Pooley, C., Sikora, L.J., 1996. Fractal parameters of pore surfaces as derived from micromorphological data: effect of long-term management practices. Geoderma 74, 305–319. Papadopoulos, A., Bird, N., Whitmore, A., Mooney, S., 2009. Investigating the effects of organic and conventional management on soil aggregate stability using X-ray computed tomography. Eur. J. Soil Sci. 60, 360–368. Peel, M.C., Finlayson, B.L., McMahon, T.A., 2007. Updated world map of the köppen-geiger climate classification. Hydrol. Earth Syst. Sci. Discuss. 4, 439–473. Perfect, E., Kay, B., 1995. Applications of fractals in soil and tillage research: a review. Soil Tillage Res. 36, 1–20. Perret, J., Prasher, S., Kacimov, A., 2003. Mass fractal dimension of soil macropores using computed tomography: from the box-counting to the cube-counting algorithm. Eur. J. Soil Sci. 54, 569–579. Pierret, A., Capowiez, Y., Belzunces, L., Moran, C., 2002. 3D reconstruction and quantification of macropores using X-ray computed tomography and image analysis. Geoderma 106, 247–271. Posadas, A.N., Giménez, D., Quiroz, R., Protz, R., 2003. Multifractal characterization of soil pore systems. Soil Sci. Soc. Am. J. 67, 1361–1369. Rachman, A., Anderson, S., Gantzer, C., 2005. Computed-tomographic measurement of soil macroporosity parameters as affected by stiff-stemmed grass hedges. Soil Sci. Soc. Am. J. 69, 1609–1616. Rasiah, V., Alymore, L., 1998. The topology of pore structure in cracking clay soil: I. the estimation of numerical density. J. Soil Sci 39, 303–314. Rogasik, H., Crawford, J., Wendroth, O., Young, I., Joschko, M., Ritz, K., 1999. Discrimination of Soil Phases by Dual Energy X-Ray Tomography. Roisin, C.J., 2007. A multifractal approach for assessing the structural state of tilled soils. Soil Sci. Soc. Am. J. 71, 15–25. Russell, D.A., Hanson, J.D., Ott, E., 1980. Dimension of strange attractors. Phys. Rev. Lett. 45, 1175–1178.

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013

ARTICLE IN PRESS 8

I. Torre et al. / Geoderma xxx (2017) xxx–xxx

Saa, A., Gascó, G., Grau, J.B., Antón, J.M., Tarquis, A.M., 2007. Comparison of gliding box and box-counting methods in river network analysis. Nonlinear Process. Geophys. 14, 603–613. Santiago, A., Cárdenas, J., Losada, J., Benito, R., Tarquis, A., Borondo, F., 2008. Multiscaling of porous soils as heterogeneous complex networks. Nonlinear Process. Geophys. 15, 893–902. Siqueira, G., Silva, E., Montenegro, A., Vidal Vázquez, E., Paz-Ferreiro, J., 2013. Multifractal analysis of vertical profiles of soil penetration resistance at the field scale. Nonlinear Process. Geophys. 20, 529–541. Soil Survey Staff, 2010. Key to Soil Taxonomy. 11th ed., USDA-Natural Resources Conservation Service, U.S. Gov. Printing Office, Washington, USA. Taina, I., Heck, R., Elliot, T., 2008. Application of X-ray computed tomography to soil science: a literature review. Can. J. Soil Sci. 88, 1–19. Tarquis, A., Giménez, D., Saa, A., Díaz, M., Gascó, J., 2003. Scaling and Multiscaling of Soil Pore Systems Determined by Image Analysis. Scaling Methods in Soil Physics. CRC Press, Boca Raton, FL., pp. 19–34. Tarquis, A., Heck, R., Andina, D., Alvarez, A., Anton, J., 2009. Pore network complexity and thresholding of 3d soil images. Ecol. Complex. 6, 230–239. Tarquis, A., Heck, R., Grau, J., Fabregat, J., Sanchez, M., Antón, J., 2008. Influence of thresholding in mass and entropy dimension of 3-d soil images. Nonlinear Process. Geophys. 15, 881–891.

Tarquis, A., Platonov, A., Matulka, A., Grau, J., Sekula, E., Diez, M., Redondo, J., 2014. Application of multifractal analysis to the study of sar features and oil spills on the ocean surface. Nonlinear Process. Geophys. 21, 439–450. Tarquis, A.M., Mèndez, V., Grau, J.B., Antòn, J.M., Andina, D., 2007. Fractals as pre-processing tool for computational intelligence application. . Computational IntelligenceSpringer., pp. 193–212. Taylor, J., 1997. An Introduction to Error Analysis: The Study of Uncertainties in Physical Measurements, 327 pp. univ. Teng, H.T., Ewe, H.-T., Tan, S.L., 2010. Multifractal dimension and its geometrical terrain properties for classification of multi-band multi-polarized SAR image. Prog. Electromagn. Res. 104, 221–237. Vázquez, E.V., Moreno, R.G., Miranda, J., Díaz, M., Requejo, A.S., Ferreiro, J.P., Tarquis, A., 2008. Assessing Soil Surface Roughness Decay During Simulated Rainfall by Multifractal Analysis. Young, I.M., Crawford, J.W., 2004. Interactions and self-organization in the soilmicrobe complex. Science 304, 1634–1637. Zhou, H., Perfect, E., Li, B., Lu, Y., 2010. Effects of bit depth on the multifractal analysis of grayscale images. Fractals 18, 127–138. Zhou, H., Perfect, E., Lu, Y., Li, B., Peng, X., 2011. Multifractal analyses of grayscale and binary soil thin section images. Fractals 19, 299–309.

Please cite this article as: I. Torre et al., Multifractal analysis of 3D images of tillage soil, Geoderma (2017), http://dx.doi.org/10.1016/j. geoderma.2017.02.013