Multifractal analysis in soil properties: Spatial signal versus mass distribution

GEODER-12428; No of Pages 12 Geoderma xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

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

Multifractal analysis in soil properties: Spatial signal versus mass distribution M.C. Morató a,⁎, M.T. Castellanos a, N.R. Bird b, A.M. Tarquis a,b,c a b c

Dpto. Matemática Aplicada, Universidad Politécnica de Madrid (UPM), Avda. Complutense s/n. 28040 Madrid, Spain C.E.I.G.R.A.M., E.T.S.I.A.A.B., UPM, Madrid, Spain Grupo de Sistemas Complejos, UPM, Avda. Complutense s/n, 28040 Madrid, Spain

a r t i c l e

i n f o

Article history: Received 22 January 2016 Received in revised form 26 July 2016 Accepted 4 August 2016 Available online xxxx Keywords: Mass function Multifractal spectrum Generalized Structure Function Scaling exponent function

a b s t r a c t The spatial variability of soil properties is a constant expected factor that must be considered in soil studies. This variability is composed of “functional” variations and random fluctuations or noise. Multifractal formalism is suitable for variables with self-similar distributions on a spatial domain. Multifractal analysis can provide insight into the spatial variability of soil parameters. In soil science, it has been quite popular to characterize the scaling property of a variable measured along a transect as a mass distribution of a statistical measure on a length domain of the studied transect. The analysed variable is divided into a number of self-similar segments, and the partition function and mass function are estimated. Based on these estimations, the multifractal spectrum (MFS) is calculated. Another technique that can be applied focuses on the variations of a measure by analysing the absolute differences in the soil property values at different scales, such as the Generalized Structure Function (GSF) and the Universal Multifractal Model (UMM). The aim of this study was to compare both types of multifractal methods on a set of soil physical properties measured on a common 1024 m transect across arable fields at Silsoe in Bedfordshire, East-Central England. The studied properties were total porosity (Porosity), gravimetric water content (GWC) and nitrous oxide flux (N2O flux). The results showed that when using both methods, the N2O flux exhibits a distinctive multifractal character, and weak multifractal characters are detected in the GWC and Porosity cases. Additionally, several parameters were calculated and discussed. Finally, the relationship between the mass exponent function (τ(q)) and the GSF (ζ(q)) found in the literature, was positively verified for the three variables. On the contrary, the relationship between ζ(q) and the scaling exponent function based on UMM (K(q)) showed discrepancies in N2O flux and GWC for q values higher than 3. This is the first time that these comparisons have been made on soil property data. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Soil properties, such as pH, soil moisture, and porosity, vary spatially and exhibit strong fluctuations even over short distances. This variability is due to the combined action of physical, chemical and biological processes that operate with different intensities and at different scales. The description and quantification of the spatial variability of soil properties are important for modelling soil processes (Burrough et al., 1994). This variability is composed of “functional” (defined) variations and random fluctuations or noise (Goovaerts, 1997, 1998). However, the distinction between these two components is scale dependent because increasing the scale of observation almost always reveals structure in the noise (Logsdon et al., 2008). Geostatistical methods and, more

⁎ Corresponding author. E-mail address: [email protected] (M.C. Morató).

recently, multifractal/wavelet techniques have been used to characterize the scaling and heterogeneity of soil properties along with other methods originating from complexity science (de Bartolo et al., 2011). Many fractal/multifractal methods have been developed to characterize these features over the years. Halsey et al. (1986) formulated the fixed-size box-counting algorithm, which is the most common classical multifractal analysis (MFA) method, to calculate the multifractal exponents, such as the scaling exponent (τ(q)) and the generalized fractal dimension D(q). This method has been widely used in many soil science studies (Folorunso et al., 1994; Kravchenko et al., 2002, 2003; Vereecken et al., 2007). Hurst (1951) proposed a rescaled range analysis (R/S analysis) to study the Nile and the problems related to water storage. More importantly, he proposed an important exponent, generally known as the Hurst exponent, to quantify the long-range correlations of the signal series. However, the R/S analysis can only handle stationary signals. To handle the fluctuations in non-stationary signals, new methods arise mainly from cascade models and turbulence studies (Davis et

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

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

2

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx 0.8

0.20

Porosity (%)

Porosity (%)

0.7 0.6 0.5 0.4 0.3 0.2

0.15 0.10 0.05

0.1 0.0

0.00 1

26

51

76

101 126 151 176 201 226 251

1

26

51

76 101 126 151 176 201 226 251

data point (i)

data point (i) 0.6 0.20

GWC (%)

GWC (%)

0.5 0.4 0.3

0.16 0.12

0.2

0.08

0.1

0.04

0.0

0.00 1

26

51

76

101 126 151 176 201 226 251

1

26

51

76 101 126 151 176 201 226 251

data point (i)

data point (i) 350 350 300

250

Ud. N2O

Ud. N2O

300

200 150 100 50

250 200 150 100 50 0

0 1

26

51

76

1

101 126 151 176 201 226 251

26 51 76 101 126 151 176 201 226 251

data point (i)

data point (i)

Fig. 1. Original data of the soil variables: Porosity (%), Gravimetric Water Content (GWC) (%) and N2O flux on the left column. On the right side, the absolute differences obtained with lag 1 of the corresponding variable.

al., 1994; Schmitt et al., 1995; Taqqu et al., 1995). Applying these methods, researchers determined the fractal scaling properties and the long-range correlations in both stationary and non-stationary series. In many soil studies, researchers have characterized the scaling property of a variable measured along a transect as a mass distribution on a spatial domain of the studied field (Zeleke and Si, 2004, 2006). For this characterization, the transect is divided into a number of self-similar segments. The differences among the subsets are identified using D(q) and a multifractal spectrum (Folorunso et al., 1994; Caniego et al., 2005; Tarquis et al., 2008a). Recently, several authors (Siqueira et al., 2013; Lopez de Herrera et al., 2016) have applied Multifractal Analysis to profiles of soil penetrometer resistance data sets and found that these methods added complementary information to describe the spatial arrangement to methods of classical statistics. However, only recent works on agricultural soils have studied the application of these methods to cases in which a measure along a transect is observed as a random signal. Pozdnyakova et al. (2005) evaluated the spatial variability of cranberry yield by applying a Generalized Structure Function, proving the influence of multiscale factors (nonlinear structure functions). Kravchenko (2008) approached the spatial features of environmental and agronomic variables using multifractal characteristics in a stochastic simulation. Garcia Moreno et al. (2010) assessed the variability of soil surface roughness using the Generalized Structure Function of transects to compare soil types and tillage tools, with promising results. In this work, we focused on the use of MFA to study the relation of the characterization of the measure among different scales. Comparing different methods, we found that there are several works studying

wider scaling behaviours which cannot be captured in a consistent way by the MFA. These analyses include extended power-law scaling (linear relations between log structure functions of successive orders) at all lags, and frequency distributions of the variables' increments, which tend to be symmetric with peaks that grow sharper and tails that become heavier as the lags between pairs of values decrease (Guadagnini et al., 2015, 2014; Riva et al., 2015). Based on the foregoing, the present study aimed to apply both types of MFA methods, the Generalized Structure Function and Multifractal Spectrum, to data on soil properties along a transect of arable fields, to compare and evaluate the results obtained for characterizing their structure and variability. 2. Material and methods 2.1. Experimental site The data used in this paper were collected in a survey on a common 1024 m transect across arable fields at Silsoe in Bedfordshire, Table 1 Descriptive statistics using the first fourth moments (average, variance, asymmetry and kurtosis) of soil porosity (Porosity), gravimetric water content (GWC) and N2O flux values (N2O). Statistics

Porosity

GWC

N2O

Average Variance Asymmetry Kurtosis

0.5736 0.0040 −0.8559 0.9440

0.3475 0.0055 −0.4289 −0.8398

54.42 2970.74 1.59 2.81

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

3

Table 2 Descriptive statistics using the first fourth moments (average, variance, kurtosis and asymmetry) of the differences in value of soil porosity (Porosity), gravimetric water content (GWC) and N2O flux values (N2O) at different lags. Statistics

Porosity (%)

Lag Average Variance Kurtosis Assimetry Data points Statistics Lag Average Variance Kurtosis Assimetry Data points Statistics Lag Average Variance Kurtosis Assimetry Data points

1 −0.0009 0.0017 0.1756 0.1317 256

2 −0.0019 0.0021 −0.0304 0.3275 128

4 −0.0037 0.0024 −0.0144 −0.4688 64

1 −0.0007 0.0007 1.5402 0.4857 256

2 −0.0014 0.001 2.4441 0.371 128

4 −0.0029 0.0014 1.8047 0.1189 64

1 0.0938 2314.4225 8.2063 0.9069 256

2 0.1875 3698.5315 7.4699 0.8156 128

4 0.375 3521.7937 1.3533 −0.6926 64

8 −0.0075 0.0038 −0.5959 −0.143 32 GWC (%) 8 −0.0058 0.0027 2.792 −1.068 32 N2O(%) 8 0.75 1887.6129 0.7432 −0.0886 32

16 −0.0149 0.0086 −0.8971 0.487 16

32 −0.0299 0.0079 0.7306 −0.6786 8

64 −0.028 0.0046 −5.6295 −0.0587 4

16 −0.0115 0.0041 0.1414 −0.5364 16

32 −0.0231 0.0075 1.7053 0.951 8

64 −0.0018 0.0047 2.7195 1.4518 4

16 1.5 3261.4667 4.0086 1.7007 16

32 3 4445.4286 3.171 1.3594 8

64 7 5776.6667 0.2131 −0.7422 4

a

L))

0.0

ln(Mq(

5.0

10.0 15.0

5.0

4.5

4.0

3.5

3.0

2.5

ln(

2.0

1.5

1.0

0a.5

0.0

L)

b

ln(Mq( i/L)

0.0 5.0 10.0 15.0 5.0

4.5

4.0

3.5

3.0

2.5

2.0

1.5

1.0

0.5

0.0

ln( i/L)

c

lnMq(( i))

20 15 10 5 0 5.0

4.5

4.0

3.5

3.0

ln(

2.5

2.0

1.5

1.0

0.5

0.0

L)

Fig. 2. Bilog plot of Generalized Structure Function (Mq(Δi)) versus ΔiΔimax being L=Δimax for: a) Porosity and b) Gravimetric Water Content (GWC) and c) N2O flux. The arrows point out the range of scale selected. The different symbols correspond to the different values of the exponent q (from top q = 0.25 to 4).

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

4

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

The same calculations were performed on each variable after differentiating the series at several non-overlapping lags, from 1 till 64. In this way, we could study the statistical moments of the frequency distribution of the values obtained in each lag. When a differentiation with lag 1 was applied, the statistics were calculated on 256 values for each variable. In the case of differentiation with lag 64, the number of values obtained was 4. For all these calculations, XLStat-Pro software (Addinsoft, New York) was used. This is a statistical package for Microsoft Excel that has been in development since 1993.

east-central England. The data has previously been described by Lark et al. (2004). The first sample point on the transect was at UK Ordnance Survey (OS) co-ordinates 508570, 235605, and the soil was sampled at 256 locations at 4 m intervals on a line running on a bearing of 188 degrees relative to UK OS grid north. The transect crossed four arable fields: Cashmore field (positions 1–50), Banqueting field (positions 57–76), Obelisk East field (positions 80–152) and Bypass field (positions 154–256). Intervening positions were in the waste ground under the grass. Cashmore field had been drilled with a winter barley crop 14 days before the soil was sampled, and Bypass and Obelisk East fields with winter wheat 12 days before. Banqueting field was in wheat stubble. This variability in the solid and drift geology influences the properties of the surface soil. The data selected from this survey for analysis in this paper were porosity, gravimetric water content and nitrous oxide flux. The values of all these variables are shown in Fig. 1 (left side).

2.3. Generalize Structure Function (GSF) and Universal Multifractal Model (UMM) The structure function analysis basically consists of studying the scaling behaviour of the non-overlapping fluctuations of a variable for different scale increments. The statistical moments of these fluctuations are estimated, which depend only on the scale increment (Monin and Yaglom, 1975). For non stationary processes GSF of order q is defined as the qth moment of the increments of initial values μ(i). The equation is:

2.2. Descriptive statistics


 Mq ðΔiÞ ≡ jμ ði þ ΔiÞ−μ ðiÞjq

For each of the variables of this study, the first four statistical moments—average, variance, kurtosis and asymmetry (skewness)—were calculated, to study their similitude with a Gaussian distribution.

ð1Þ

where i denotes the ith data point, and 〈〉 denotes the ensemble average. GSF are generalized correlation functions, which is particularly evident

a 2.5 2.0 1.5 1.0 0.5 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

q

b 2.5 2.0 1.5 1.0 0.5 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

2.5

3.0

3.5

4.0

4.5

q

c

2.5 2.0 1.5 1.0 0.5 0.0 0.0

0.5

1.0

1.5

2.0

q Fig. 3. ζ(q) curves obtained by Generalized Structure Function (GSF) for: a) Porosity, b) Gravimetric Water Content (GWC) and c) N2O flux. Continuous lines are straight lines with slope 0.5 representing a non correlated noise.

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

a

5

1.0

H(q)

0.8 0.6 0.4 0.2 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

2.5

3.0

3.5

4.0

4.5

q

b 1.0

H(q)

0.8 0.6 0.4 0.2 0.0 0.0

0.5

1.0

1.5

2.0

q

c

1.0

H(q)

0.8 0.6 0.4 0.2 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

q Fig. 4. Generalized Hurst exponent (H(q)) for: a) Porosity, b) Gravimetric Water Content and c) N2O flux. Continuous lines are straight lines with slope 0 corresponding to non correlated noise with Hurst value 0.5.

from Eq. (1) for the case of q = 2. In general, q may be any real number not just integers, and can even be negative. However, there are divergence problems inherent to the negative-order exponent so that computations are best restricted to positive real number (Davis et al., 1994). If the process μ(i) is scale-invariant and self-similar or self-affine over some range of space lags Δimin ≤ Δi ≤ Δimax, then the qth-order structure function is expected to scale as: ζ ðqÞ

Mq ðΔiÞ ¼ C q Δi

≈Δi

ζ ðqÞ

ð2Þ

where Cq can be a function of Δi which varies more slowly than any power of Δi, and ζ(q) is the exponent of the structure function. ζ(q) is a monotonically non-decreasing function of q if μ(i) has absolute bounds (Frisch, 1995; Marshak et al., 1994). From Eq. (2) we can see that the statistics of the fluctuations over space lags Δi has two components; the first because it depends of the fluctuations at low Δivalues, the second because of the scaling relation between the fluctuations and Δi. From ζ(q) the first moment is related to the degree of non conservation of a given field, H = ζ(q =1).

The behaviour described by Eqs. (1) and (2) is called “multiscaling” because each statistical moment is scaling but with a different exponent. Therefore, a hierarchy of exponents can be defined using ζ(q): H ðqÞ ¼

ζ ðqÞ q

ð3Þ

where H(q) is the generalized Hurst exponent or self-similarity scaling exponent (Davis et al., 1994). Calculation of H(q) allows the

Table 3 Generalized Hurst exponents values (H(q)) from the two extremes, q = 0.25 and q = 4, and for q = 1 (ζ(1) = H(1) = H) for Porosity, Gravimetric Water Content (GWC) and N2O flux. ΔH is the amplitude of variation of H(q)(ΔH=H(0.25)−H(4)). Estimations have been made through the generalized structure function (GSF).

Porosity (%) GWC (%) N2O flux (Ud.)

H(0.25)

H(4)

ΔH

H(1)

0.3611 0.4466 0.5124

0.225 0.3235 0.3264

0.1360 0.1231 0.1860

0.3117 0.4072 0.4516

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

6

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

straightforward identification of persistence, or long-space correlation, as well as the stationary/non stationary and monofractal/multifractal nature of the data (Lovejoy et al., 2001). Stationary processes have scale-independent increments and ζ(q) = H(q) ≡ 0, due to the invariance under translation. Processes with a linear ζ(q) (or a constant H(q)) are non-stationary and monofractal, otherwise they are non-stationary and multifractal. The UMM assumes that multifractals are generated from a random variable with an exponentiated extreme Levy distribution (Lavallée et al., 1991; Tessier et al., 1993). In UM analysis, the scaling exponent K(q) is highly relevant. K(q) relates the moment of order q of the field with the scale by the next expression:  K ðqÞ Δ Mq ðΔi Þ ¼ ðΔμ ðΔi ÞÞ ≈ i max Δi


q

Dq ¼

1 logχ ðq; δÞ lim q−1 δ→0 logδ

ð9Þ

nðδÞ

where χ(q, δ) is the partition function (PF) defined as χ(q, δ) = ∑ μ qi i¼1

ðδÞ, q is a real number, −∞ b q b ∞, and n(δ) the number of intervals of size δ (Chhabra and Jensen, 1989). A log-log plot of the quantityχ(q, δ) versus δ for different values of q yields: χ ðq; δÞ∞δτðqÞ

ð10Þ

ð4Þ

where ΔiΔimax is the scale ratio, which is inversely proportional to the size of the measurement interval. The scaling exponent function K(q) for the moments q of a cascade conserved process is obtained according to Schertzer and Lovejoy (1987) as follows: 8 9 < C 1 ðqαL −qÞ = if α L ≠1 K ðqÞ ¼ α −1 L : ; C 1 q logðqÞ if α L ≡ 1

defined by Hentschel and Procaccia (1983) as:

ð5Þ

where τ(q) is the mass exponent or correlation exponent of the qth order (Halsey et al., 1986; Vicsek, 1992) defined as: τðqÞ ¼ ðq−1ÞDq

a

ð11Þ

0.8 0.7 0.6 0.5

k(q)

where C1 is the mean intermittency codimension and αL is the Levy index. These are known as the UM parameters. To obtain these multifractality parameters, the Double Trace Method (DTM) was applied. For further details see Lavallée et al. (1991) and Tessier et al. (1993).

0.4 0.3 0.2 0.1 0.0

2.4. Generalized dimensions and multifractal Spectrum

0.1 0.0

Multifractal models provide more information about a distribution of a physical system than fractal models (Voss, 1988). The probability mass function, μi(δ) describes the portion of total mass contained in each segment of the transect, and was estimated as:

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

q

b

0.8 0.7 0.6

ð6Þ

where Mi(δ) is the value of the measure in the i-th segment of size δ, and Mt. represents the total mass in the whole transect. For heterogeneous or non uniform systems the probability mass function, μi(δ), scales with the interval size, δ, as:

0.5

k(q)

μ i ðδÞ ¼ M i ðδÞ=Mt

0.4 0.3 0.2 0.1 0.0

MðaÞ∝δ

0.5

1

1.5

2

2.5

3

3.5

4

4.5

q

where αi is the singularity or Lipschitz–Hölder or Hölder exponent characterizing scaling in the ith region. The Hölder exponent may be interpreted as a crowding index for the degree of concentration of the measure, μ. Also, for multifractal distributed measures, the number Mδ(α) of segments of size δ, having a singularity or Hölder exponent equal to α, obeys a power law: − f ðα Þ

0.1 0

ð7Þ

ð8Þ

c 0.8 0.7 0.6 0.5

k(q)

μ i ðδÞ∝δαi

0.4 0.3 0.2 0.1

where the exponent f(α) can be defined as the fractal dimension (singularity spectrum). It describes the statistical distribution of the singularity exponent α, or in other words, counts how often specific values α of the singularity strengths occur (Feder, 1988). Multifractal sets can also be characterized on the basis of the generalized dimensions, Dq, of the qth moment orders of a distribution,

0.0 0.1 0.0 0.5 1.0

1.5 2.0 2.5 3.0

3.5 4.0 4.5

q Fig. 5. K(q) curves obtained by the Universal Multifractal Model (UMM) for: a) Porosity and b) Gravimetric Water Content (GWC), c) N2O flux.

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

The connection between the power exponents f(α) and τ(q) is made through a Legendre transformation: f ðα ðqÞÞ ¼ qα ðqÞ−τðqÞ and α ðqÞ ¼

N ðδÞ X

f ðα ðqÞÞ∝

dτðqÞ dq

2.5. Relation between ζ(q), K(q) and τ(q)

ð13Þ

1

There are several methods for multifractal scaling pattern characterization, besides the two we are describing here. A quite complete review of these can be found in the book of Seuront (2010), where detailed explanations of UMM and GSF are given. Both methods were initially developed in the area of turbulence. The space variability of a measure along a transect can often be broken up into various “scaling ranges” over which the fluctuations vary in

Singularity strength α(q) and singularity spectrum f(α):

α ðqÞ∝ i¼1

μ i ðq; δÞ log½μ i ðδÞ ð14Þ

logðδÞ

ð15Þ

logðδÞ

The graph of f(α) versus α called the multifractal spectrum, typically has a parabolic concave downward shape, with the range of α values increasing with the increase in the heterogeneity of the measure. An homogeneous fractal exhibits a narrow f(α) spectrum.

There are several ways to calculate the f(α) singularity spectrum. Following the general methodology of Evertsz and Mandelbrot (1992) and Feder (1988) and the specific techniques of Chhabra and Jensen (1989) the normalized measure μi(q,δ) that correspond to contributions of individual segments of the partition function is defined by

N ðδÞ X

μ i ðq; δÞ log½μ i ðq; δÞ

i¼1

ð12Þ

nðδÞ X μ qi ðδÞ μ i ðq; δÞ ¼ μ i q ðδÞ=

7

ln

a 35 30 25 20 15 10 5 0 5 0 10 15 20

1

2

3

4

5

6

ln

b 35 30 25 20 15 10 5 0 5 0 10 15 20

1

2

3

4

5

6

ln

c 55 50 45 40 35 30 25 20 15 10 5 0 5 0 10 15 20

1

2

3

4

5

6

Fig. 6. Bilog plot of Partition Function (χ(q,δ)) versus number of intervals size (δ) for: a) Porosity, b) Gravimetric Water Content and c) N2O flux. The different symbols correspond to the different values of the exponent q (from top q = 0.25 to 4).

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

8

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

a power law manner with respect to scale. Over these ranges, the fluctuations follow   Δμ ¼ φΔi ΔiH ;

ð16Þ

where Δi denotes the lag of the chosen soil transect, φ denotes the turbulent flux and H is the degree of non-conservation mentioned earlier. Taking the qth power of Eq. (16) and averaging over the ensemble, we see that the statistical characterization of the fluctuations is Lovejoy and Schertzer (2012)

 Mq ðΔiÞ ¼ ðΔμ ðΔiÞÞq ¼ φqΔi ΔiqH ¼

  Δi max K ðqÞ qH Δi ≈Δiζ ðqÞ ; Δi

In the case that the field is intermittent – for example, if it is the result of a multifractal process – then the exponent K(q) is generally non-linear, convex and characterizes the intermittency. Because the ensemble mean of the flux spatially averaged at any scale is the same (i.e., 〈ϕΔi〉 is a constant independent of Δi), we have K(1) = 0 and ζ(1) = H. In addition, the relationship between the UMM and the multifractal formalism based on τ(q) is described by the follow equation (Gagnon et al., 2003; Aguado et al., 2014): τðqÞ ¼ ðq−1ÞE−K ðqÞ:

ð19Þ

E is the Euclidean dimension in which the measure is embedded. We can substitute K(q) using Eq. (18), obtaining:

ð17Þ

τðqÞ ¼ ðq−1ÞE−qH þ ζ ðqÞ:

ð20Þ

where K(q) is the scaling moment function that characterizes the flux (Schertzer and Lovejoy, 1985) and is explained earlier in the UMM. Eq. (17) relates ζ(q) with K(q), and the latter could be interpreted as an intermittent correction; hence, it expresses the deviation of ζ(q) from linearity due to intermittency (Seuront et al., 1999):

In our case, E = 1 as we are working with soil transect data, giving the expression:

ζ ðqÞ ¼ qH−K ðqÞ or K ðqÞ ¼ qH−ζ ðqÞ or K ðqÞ þ ζ ðqÞ ¼ qH:

The relation between ζ(q) and K(q) (Eq. (18)) has been discussed in a recent paper (Renosh et al., 2015). It is shown that such a

a

ð18Þ

τðqÞ ¼ ðq þ ζ ðqÞÞ−ð1 þ qHÞ:

ð21Þ

6

(q)

3 0 -3 -6 -9 -6

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

q

b

6

(q)

3 0 -3 -6 -9 -6

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

q

c

6

(q)

3 0 -3 -6 -9 -6

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

q Fig. 7. Mass function (τ(q)) curves obtained based on Partition Function for: a) Porosity, b) Gravimetric Water Content (GWC) and c) N2O flux.

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

a f( )

1.0 0.8 0.6 0.4 0.6

b f( )

relation is valid for lower order moments, but for q N 2 it is not valid in general. The authors note that the problem comes from the fact that ϕΔi in Eq. (16) is defined as a small-scale positive quantity, obtained from the small-scale gradient of the field μ. However, the upscaling of ϕΔi to a given scale is not necessarily directly related to the gradient of the original field μ at the same scale. In other words, the large-scale average of an absolute value is not the same as the absolute value at large scale. The relation [21] will be tested with the results obtained from the mass exponent, based on the Partition Function, and the GSF. With respect to the relation between GSF and UMM, Eq. (18) will be analysed by plotting ζ(q) + K(q)and qH versus q, following the test applied by Renosh et al. (2015).

9

0.8

1.0

1.2

1.4

1.6

1.0 0.8 0.6

3. Results and discussion 0.4

3.1. Descriptive statistics

3.2. Generalized structure function and universal multifractal model As reported in other works (Tarquis et al., 2008b), at higher q values the errors in the GSF are higher. For this reason, we performed the analysis for a range of q values from 0.5 till 4 with increments of 0.25. The GSFs obtained from the three variables studied are shown in Fig. 2, where the maximum increment chosen was 128 data points (L = Δimax), which is equivalent to 512 m. We can observe that the bilog plot of Mq(Δi)does not show linear behaviour in every interval. For this reason, a range of ΔiΔimax intervals was selected for each variable, as in-

dicated by the arrows, obtaining a minimum R2 of 0.97 for all q values used. The results obtained from the ζ(q) curves, shown in Fig. 3, corroborate the multifractal nature of N2O flux, and much more weakly corroborate the multifractal nature of Porosity, which looks almost like a straight line in this plot. GWC exhibits a more relevant curvature in ζ(q), compared to Porosity. For comparison, we also plotted the

0.8

1.0

1.2

1.4

1.6

c 1.0

f( )

The first four statistical moments were calculated for the three variables (Table 1). From the original values, the variable with the highest variance was N2 O flux, as expected after observing their values in Fig. 1 (left side, showing the original series of μ i ) to be much lower in the Porosity and GWC cases. With respect to higherorder moments, either asymmetry or kurtosis are closer to values corresponding to a normal distribution. In the case of Porosity and N2O flux, the kurtosis is positive, indicating a “peaked” distribution. Meanwhile, GWC data showed a negative kurtosis, indicating a “flat” distribution. Studying the measured asymmetry, Porosity and GWC were both negative and therefore were slightly skewed left; meanwhile, N2 O flux had a higher positive value, indicating that the distribution is skewed right. The differentiation with lag 1 of the absolute values are visualize for each variable in Fig. 1 (right side, showing the | μi − μi + 1 | series). Focusing our attention of the differentiate values at different lags (Table 2), we can observed that lags bigger than 8 don't have enough data points to have a good estimation in the statistical moments of the frequencies distribution and therefore we will concentrated in lags from 1 to 8. In the three variables all the averages have a value close to zero, presenting Porosity and GWC a negative sign and N2O positive one. As the lag increases, Porosity and GWC decrease the average value meanwhile N2O increases. Respect to kurtosis, N2O presents values higher than Gaussian distribution that decrease as the lag increases. GWC presents kurtosis values lower than N2O and the same tendency with the lag meanwhile Porosity has values very close to a Gaussian distribution. In the case of asymmetry, all the variables and lags from 1 to 8 show values close to zero except for GWC at lag 8 that presents an asymmetry close to 1.

0.6

0.8 0.6 0.4 0.6

0.8

1.0

1.2

1.4

1.6

Fig. 8. Multifractal spectrum (f(α) versus α) for: a) Porosity, b) Gravimetric Water Content (GWC) and c) N2O flux.

line ζ ðqÞ ¼ 2q , which corresponds to Brownian motion. In the three cases, the ζ(q) function is significantly different from the line corresponding to Brownian motion, demonstrating the correlation of the increments among several scales. The generalized Hurst exponents derived from the ζ(q) function based on GSF are shown in Fig. 4, where the straight line represents the case of pure noise or uncorrelated noise (H(q) = 0.5). At lower q values the H(q) is closer to 0.5; however, H(q) decreases in the higher moments, attaining an anti-persistent character. To quantify the variation of H(q) versus q, we have subtracted the extreme values (H(0.25)− H(4)) for each variable (see Table 3). The highest amplitude of H(q) corresponds to N2O flux, which showed a more pronounced curvature in ζ(q) (see Fig. 3). The three studied series exhibit degrees of non-conservation that are lower than one (H =ζ(q = 1)=H(q = 1)) and therefore it is not necessary to apply a detrending fluctuation analysis (Kantelhardt et al., 2002). At this point, we would like to remark that the ζ(q) function is based on the behaviour of the increments of the variables at different lags. However, these increments are in absolute value, implying a lack of ability to interpret the statistics of the probability density function of the increments as other scaling models do (Riva et al., 2015).

Table 4 Parameters extracted from the multifractal spectrum (MFS) for the three variables studied: Porosity, Gravimetric Water Content (GWC) and N2O flux. Holder exponent at q = 1 (α(1)), at q = +4 (αmin), at q = −4 (αmax) and αmax − αmin (Δα). Hausdorff dimension at αmin(f(αmin)) at αmax(f(αmax)) and f(αmax) – f(αmin) (Δf).

Porosity (%) GWC (%) N2O flux (Ud.)

α(1)

αmin

αmax

Δα

f(αmin)

f(αmax)

Δf

0.9984 0.9938 0.9740

0.9922 0.9768 0.7892

1.0248 1.0734 1.8812

0.0326 0.0966 1.0921

0.9847 0.9584 0.6586

0.9462 0.8682 0.0532

−0.0385 −0.0902 −0.6053

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

10

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

All the MFS were concave down parabolic curves (Fig. 8) with a variable symmetry depending on the studied variable. We recall that in a fractal system, the MFS was shown as a single point; therefore, our results support the hypothesis of multifractal behaviour. However, in the cases of Porosity and GWC, the MFS is much more constrained than in the N2O flux case, which is a consequence of the linearity of the τ(q) function. Now we focus on the quantification of MFS amplitudes and symmetries through the difference of the extreme singularities (Δα = αmax − αmin) and the difference of their respective f(α) values (Δf = f(αmax) − f(αmin), as shown in Table 4 for the three variables. In this way, higher Δα indicates higher complexity of the structure studied in the transect; higher Δf indicates higher asymmetry in the MFS (left-handed if Δf = 0 and right-handed if Δf N 0). Studying the evolution of Δα among the three variables, we can see that there is an increase from Porosity to N2O flux, indicating an increasing complexity in the hierarchical spatial structure among scales. With respect to the behaviour of Δf in each variable (Table 4), it is again N2O flux that exhibits a higher asymmetry, showing the MFS a slight

The results obtained from the K(q) curves (Fig. 5) indicate, in agreement with the results obtained from the ζ(q) curves, the multifractal character of the three variables studied. N2O flux showed the sharpest increase in K(q). The increase was much lower for Porosity. These K(q) curves will be used to test Eq. (18). 3.3. Multifractal Spectrum A multifractal analysis was performed for each of the series obtained. In all cases we found a linear relationship between the log-log plot of the partition function (χ(q, δ)) versus the length of the scale (δ), obtaining R2 higher than 0.98 (see Fig. 6). From the estimated slope for each mass exponent (q), a τ(q) function is obtained (Fig. 7). For q = 1, we obtained τ(q = 1) = 0, indicating the conservative character of the variable. Only in the case of N2O flux is this function clearly nonlinear, reflecting a hierarchical structure from one scale to the other. Therefore, we estimated the multifractal spectrum (MFS) in an interval of q = ±4 with increments of 0.5 (Fig. 8).

a

y =1.1058x + 0.008 R2 = 0.9995

4.0

(q) [MFS]

3.0 2.0 1.0 0.0 1.0 1.0

0.5

0.0

0.5

1.0

1.5

2.0

2.5

3.0

(q)[GSF] y =1.0864x + 0.0034

b

R2 = 0.9994

(q) [MFS]

4.0 3.0 2.0 1.0 0.0 1.0 1.0

0.5

0.0

0.5

1.0

1.5

2.0

2.5

3.0

(q)[GSF]

(q [MFS])

c

y = 0.9957x - 0.0055 R2 = 1

4.0 3.0 2.0 1.0 0.0 1.0 1.0

0.5

0.0

0.5

1.0

1.5

2.0

2.5

3.0

(q)[GSF] Fig. 9. Relation between estimated τ(q) through Generalized Structure Function (GSF) following Eq. (19) and τ(q) based on the Partition function (Eq. (9)) used in the Multifractal Spectrum (MFS) for: a) Porosity, b) Gravimetric Water Content (GWC) and c) N2O flux.

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

left-handed shape. Porosity and GWC are almost symmetrical. Finally, Porosity and GWC have an entropy dimension α(1) higher than 0.99, while N2O flux has a value of 0.97.

11

those of Renosh et al. (2015), where in their case the discrepancies were found for q N 2. 4. Conclusions

3.4. Relation between ζ(q), K(q) and τ(q) The relationship expressed in Eq. (21) has been tested by plotting both sides of the equation. For this purpose, we have used q values from 0.25 till 4 with an increment of 0.25, to have enough data points to check the correlation. Using the partition function, all the scale ranges were used. For a shorter range of scales, GSF showed linear behaviour in the bilog plots used to estimate the parameters. GWC exhibits a linear pattern from increments of 4 till 256 m. Porosity showed this behaviour from 4 till 64 m. In case of N2O flux, this behaviour was exhibited only from 32 to 256 m. Results are shown in Fig. 9, where for the three variables a straight line is obtained with correlations higher than 0.99. This is remarkable considering that the available data set was based on 256 points. We were expecting an unbiased estimation of the multifractal parameters (Roberts and Cronin, 1996). The relationship expressed in Eq. (18) has been tested by plotting two functions, K(q) + ζ(q) and qH, versus q. If these two functions are identical then Eq. (18) is verified with these experimental data (see Fig. 10). In the case of Porosity, for all the q tested, the relationship holds because both functions are a common straight line. When we observed the behaviour in GWC, in both functions a tendency to mutually divert exists for q larger than 3, while in N2O flux this tendency is more clear. The stronger the multifractal character of the variable, the greater the discrepancies found for q N 3. These results are in agreement with

(q)+k(q)

a

2.0 1.5 1.0 0.5

The scaling properties of the three variables (total porosity, gravimetric water content and nitrous oxide flux) along a soil transect were studied, using two techniques. Applying general structure function analysis, ζ(q) curves showed the multifractal nature of N2O flux and much more weakly the multifractal natures of Porosity and GWC. A Hurst exponent H(q) lower than 0.5 indicates the anti-persistent character of the data, mainly in the higher moments. The K(q) curves agree with the previous analysis of the multifractal characters of the three variables. Multifractal spectrum showed concave-down parabolic curves, which supports the hypothesis of multifractal behaviour. The highest multifractality corresponded to N2O flux. In addition, N2O flux had an entropy dimension of 0.97, while Porosity and GWC had values very close to 1.0. Finally, for Porosity, Gravimetric Water Content (GWC) and N2O flux, straight line relations were obtained between the τ(q) estimated through Generalized Structure Function (GSF) and the τ(q) based on the Partition function used in the Multifractal Spectrum (MFS), in agreement with theoretical studies. In contrast, the relation between the ζ(q) and K(q) curves was only found for Porosity, which is the variable with the lowest multifractal character. The relation was not obtained for GWC and N2O for moments (q) higher than 3. There are several multifractal methods that can be applied to the characterization of scaling properties. Only two types have been discussed here, which provided similar results in this study on the multifractal characters of the three variables. However, the relations found in the literature among τ(q), ζ(q) and K(q) are not always found in the data sets for several reasons. Many times, in soil studies, we have a finite data set (Roberts and Cronin, 1996), soils are not ideal multifractals (Kravchenko and Pachepsky, 2004), and there are always various uncertainties in defining the linear ranges (Tarquis et al., 2008c).

0.0 0

1

2

3

4

5

Acknowledgements

q

(q)+k(q)

b

The authors would like to thank Dr. Murray Lark of BSC for useful comments and discussions and this data set donated for this study. Our gratitude to the anonymous reviewers and editor, their suggestions has improved this work. This research has been partially supported by funding from MINECO under contract No. MTM2015-63914-P and CICYT PCIN-2014-080.

2.0 1.5 1.0 0.5

References

0.0

0

1

2

3

4

5

3

4

5

q

(q)+k(q)

c

2.0 1.5 1.0 0.5 0.0 0

1

2

q Fig. 10. Relation between the addition of ζ(q) based on the Generalized Structure Function and K(q) based on the Universal Multifractal Model (UMM) versus q: a) Porosity, b) Gravimetric Water Content (GWC) and c) N2O flux. Continuous straight line represents qH versus q.

Aguado, P.L., Del Monte, J.P., Moratiel, R., Tarquis, A.M., 2014. Spatial characterization of landscapes through multifractal analysis of DEM. Sci. World J. 2014, 563038. http:// dx.doi.org/10.1155/2014/563038 9 pp. Bartolo, S.d., Otten, W., Cheng, Q., Tarquis, A.M., 2011. Modeling soil system: complexity under your feet. Biogeosciences 8, 3139–3142. Burrough, P.A., Bouma, J., Yates, S.R., 1994. The state of the art in pedometrics. Geoderma 62, 311–326. Caniego, F.J., Espejo, R., Martín, M.A., San José, F., 2005. Multifractal scaling of soil spatial variability. Ecol. Model. 182, 291–303. Chhabra, A.B., Jensen, R.V., 1989. Direct determination of the ƒ(α) singularity spectrum. Phys. Rev. Lett. 62, 1327. Davis, A., Marshak, A., Wiscombe, W., Cahalan, R., 1994. Multifractal characterizations of nonstationarity and intermittency in geophysical fields: observed, retrieved, or simulated. J. Geophys. Res. 99, 8055–8072. Evertsz, C.J.G., Mandelbrot, B.B., 1992. In: Peitgen, H.O., Jurgens, H., Saupe, D. (Eds.), Chaos and Fractals: New Frontiers of Science. Springer-Verlag, New York, pp. 921–953. Feder, J., 1988. Fractals. Plenum Press, New York 283 pp. Folorunso, O.A., Puente, C.E., Rolston, D.E., Pinzón, J.E., 1994. Statistical and fractal evaluation of the spatial characteristics of soil surface strength. Soil Sci. Soc. Am. J. 58, 284–294. Frisch, U., 1995. Turbulence. The Legacy of a. Kolmogorov. Cambridge University Press, Cambridge, UK.

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004

12

M.C. Morató et al. / Geoderma xxx (2016) xxx–xxx

Gagnon, J.S., Lovejoy, S., Schertzer, D., 2003. Multifractal surfaces and terrestrial topography. Europhys. Lett. 62, 801–807. Garcia Moreno, R., Diaz Alvarez, M.C., Saa Requejo, A., Valencia Delfa, J.L., Tarquis, A.M., 2010. Multiscaling analysis of soil roughness variability. Geoderma 160 (1), 22–30 30. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford Univ. Press, New York. Goovaerts, P., 1998. Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties. Biol. Fertil. Soils 27, 315–334. Guadagnini, A., Neuman, S.P., Schaap, M.G., Riva, M., 2014. Anisotropic statistical scaling of soil and sediment texture in a stratified deep vadose zone near Maricopa, Arizona. Geoderma 214-215, 217–227. http://dx.doi.org/10.1016/j.geoderma.2013.09.008. Guadagnini, A., Neuman, S.P., Nan, T., Riva, M., Winter, C.L., 2015. Scalable statistics of correlated random variables and extremes applied to deep borehole porosities. Hydrol. Earth Syst. Sci. 19, 729–745. http://dx.doi.org/10.5194/hess-19-729-2015. Halsey, T.C., Jensen, M.H., Kadanoff, L.P., Procaccia, I., Shraiman, B.I., 1986. Fractal measures and their singularities: the characterization of strange sets. Phys. Rev. A 33, 1141–1151. Hentschel, H.G.E., Procaccia, I., 1983. The infinite number of generalized dimensions of fractals and strange attractors. Phys. D 8, 435–444. Hurst, H.E., 1951. Trans. Am. Soc. Civ. Eng. 116–770. Kantelhardt, J.W., Zschiegner, S.A., Bunde, A., Havlin, S., Koscielny-Bunde, E., Stanley, H.E., 2002. Multifractal detrended fluctuation analysis of nonstationary time series. Phys. A 316, 87–114. Kravchenko, A.N., 2008. Stochastic simulations of spatial variability based on multifractal characteristics. Vadose Zone J. 7 (2), 521–524. Kravchenko, A.N., Pachepsky, Y.A., 2004. Soil Variability Assessment with Fractal Techniques. In: Alvarez-Benedi, J., Munoz-Carpena, R. (Eds.), Soil-Water-Solute Process Characterization: An Integrated Approach. CRC Press, Boca Ratón, pp. 617–638 ISBN: 978-1-56670-657-5. Kravchenko, A.N., Omonode, R., Bollero, G.A., Bullock, D.G., 2002. Quantitative mapping of soil drainage classes using topographical data and soil electrical conductivity. Soil Sci. Soc. Am. J. 66, 235–243. Kravchenko, A.N., Thelen, K.D., Bullock, D.G., Miller, N.R., 2003. Relationship among crop grain yield, topography, and soil electrical conductivity studied with crosscorrelograms. Agron. J. 95, 1132–1139. Lark, R.M., Milne, A.E., Addiscott, T.M., Goulding, K.W.T., Webster, C.P., O'Flaherty, S., 2004. Scale and location-dependent correlation of nitrous oxide emissions with soil properties: an analysis using wavelets. Eur. J. Soil Sci. 55, 611–627. Lavallée, D., Schertzer, D., Lovejoy, S., 1991. On the Determination of the Codimension Function. In: Schertzer, D., Lovejoy, S. (Eds.), Non-Linear Variability in Geophysics. Springer Netherlands, Dordrecht, The Netherlands, pp. 99–109. Logsdon, S.D., Perfect, E., Tarquis, A.M., 2008. Multiscale soil investigations: physical concepts and mathematical techniques. Vadose Zone J. 7, 453–455. Lopez de Herrera, J., Herrero Tejedor, T., Saa-Requejo, A., Tarquis, A.M., 2016. Effects of tillage on variability in soil penetration resistance in an olive orchard. Soil Res. 54 (2), 134–143. Lovejoy, S., Schertzer, D., 2012. Haar wavelets, fluctuations and structure functions: convenient choices for geophysics. Nonlinear Process. Geophys. 19, 513–527. Lovejoy, S., Schertzer, D., Stanway, J.D., 2001. Fractal behaviour of ozone, wind and temperature in the lower stratosphere. Phys. Rev. Lett. 86, 5200–5203. Marshak, A., Davis, A., Cahalan, R.F., Wiscombe, W.J., 1994. Bounded Cascade models as non-stationary multifractals. Phys. Rev. E 49, 55–69.

Monin, A.S., Yaglom, A.M., 1975. Stratistical Fluid Mechanics: Mechanics of Turbulence. MIT Press, Boston. Pozdnyakova, L., Giménez, D., Oudemans, P.V., 2005. Spatial analysis of cranberry yield at three scales. Agron. J. 97 (1), 49–57. Renosh, P.R., Schmitt, F.G., Loisel, H., 2015. Scaling analysis of ocean surface turbulent heterogeneities from satellite remote sensing: use of 2D structure functions. PLoS One 10, e0126975. http://dx.doi.org/10.1371/journal.pone.0126975. Riva, M., Neuman, S.P., Guadagnini, A., 2015. New scaling model for variables and increments with heavy-tailed distributions. Water Resour. Res. 51, 4623–4634. http://dx. doi.org/10.1002/2015WR016998. Roberts, A.J., Cronin, A., 1996. Unbiased estimation of multi-fractal dimensions of finite data sets. Phys. A 233 (3–4), 867–878. Schertzer, D., Lovejoy, S., 1985. The Dimension and Intermittency of Atmospheric Dynamics. In: Launder, B. (Ed.)Turbulent Shear Flow vol. 4. Springer, New York, pp. 7–33. Schertzer, D., Lovejoy, S., 1987. Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes. J. Geophys. Res. 92, 9693–9714. Schmitt, F., Lovejoy, S., Schertzer, D., 1995. Multifractal analysis of the Greenland ice-core project climate data. Geophys. Res. Lett. 22 (13), 1689–1692. Seuront, L., 2010. Fractals and Multifractals in Ecology and Aquatic Science. CRC press, Boca Ratón 344 pp. Seuront, L., Schmitt, F., Lagadeuc, Y., Schertzer, D., Lovejoy, S., 1999. Universal multifractal analysis as a tool to characterize multiscale intermittent patterns: example of phytoplankton distribution in turbulent coastal waters. J. Plankton Res. 21 (5) 877-822. Siqueira, G.M., Silva, E.F.F., Motenegro, A.A.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. Taqqu, M.S., Teverovsky, V., Willinger, W., 1995. Estimators for long-range dependence: an empirical study. Fractals 3, 785. Tarquis, A.M., Castellanos, M.T., Morató, M.C., Bird, N.R., 2008a. Multifractal analysis on soil transect data. Geophys. Res. Abstr. 10 EGU 2008-A-04453. Tarquis, A.M., Morató, M.C., Castellanos, M.T., Perdigones, A., 2008b. Comparison of structure function and detrended fluctuation analysis in wind time series. Riv. Nuovo Cimento C 31 (5–6), 633–651. Tarquis, A.M., Heck, R.J., Grau, J.B., Fabregat, J., Sanchez, M.E., Anton, J.M., 2008c. Influence of thresholding in mass and entropy dimension of 3-D soil images. Nonlinear Process. Geophys. 15 (6), 881–891. Tessier, Y., Schertzer, D., Lovejoy, S., 1993. Universal multifractal theory and observation for rain and clouds. J. Appl. Meteorol. 32, 223–250. Vereecken, H., Kasteel, R., Vanderborght, J., Harter, T., 2007. Upscaling hydraulic properties and soil water flow processes in heterogeneous soils: a review. Vadose Zone J. 6, 1–28. Vicsek, T., 1992. Fractal Growth Phenomena. second ed. Word Scientific Publishing Co., Singapore. Voss, R.F., 1988. Fractals in Nature: From Characterization to Simulation. In: Peitgen, H.O., Saupe, D. (Eds.), The Science of Fractal Images. Springer-Verlag, New York, NY, pp. 21–69. Zeleke, T.B., Si, B.C., 2004. Scaling properties of topographic indices and crop yield: multifractal and joint multifractal approaches. Agron. J. 96, 1082–1090. Zeleke, T.B., Si, B.C., 2006. Characterizing scale-dependent spatial relationships between soil properties using multifractal techniques. Geoderma 134, 440–452.

Please cite this article as: Morató, M.C., et al., Multifractal analysis in soil properties: Spatial signal versus mass distribution, Geoderma (2016), http://dx.doi.org/10.1016/j.geoderma.2016.08.004