Application of Multifractal and Joint Multifractal Analysis in Examining Soil Spatial Variation: A Review

1. Introduction

Soil varies considerably from location to location [1] and the understanding of this variability has important applications in agriculture, environmental sciences, hydrology and earth sciences. For example, information about soil spatial variation is necessary for precision agriculture [2], environmental prediction [3], soil-landscape process modelling [4, 5], soil quality assessment [6, 7] and natural resources management. The quantification and characterization of soil spatial variability did not start until the latter half of the last century. Over the last four decades, systematic studies have identified the following characteristics of soil spatial variation.

• Spatial autocorrelation: Soil is a function of various environmental factors including climate, living organisms, relief, parent material, and time [8]. The individual or combined influence of these factors and various physical, chemical and biological processes produce different types of soil. Similar environmental factors and soil forming processes tend to have occurred at locations in close proximity to one another. It follows that soil properties measured at adjacent locations are likely to be more similar than properties measured at places located far apart [9, 10]. This is known as spatial autocorrelation or similarity.

• Scale dependence: The soil forming factors and processes can operate at different intensities and scales [4,11]. Soil biological processes or the activity of micro- and macro-organisms can occur at very small scales affecting the formation of soils and the variability in soil properties. In contrast, atmospheric, geologic and climatic variability can determine the formation of soil and the variability of soil properties over a large area. The scale dependence of the environmental factors determines that soil spatial variability is also scale dependent.

• Periodicity: Sometimes the pattern of variation in soil attributes may repeat at different scales or over a certain distance beyond its spatial autocorrelation. The repetition or quasi-cyclic variation in soil can arise from the repeated features in topography, geology or parent material, tillage operation or cultivation practices [12-14]. This cyclic behaviour is known as the periodicity.

• Nonstationarity: The linear or nonlinear spatial trends in the controlling factors and processes may result in a change in soil properties that is gradual and predictable in space. This is referred to as the nonstationarity, which can arise from the effects of topography, lithology, parent material, climate and vegetation [15]. Non-stationarity can also occur where there are distinct strata in the variation, such as different types of soil. The mean and/or variance of soil properties in one stratum may differ from that of another.

• Nonlinearity: Often the effect from different factors and processes are non-additive in nature [16] and do not follow the principle of superposition. These are the characteristics of a nonlinear system [17], which can be explained by a simple example. Soil water storage is controlled by a number of factors [such as elevation, texture, vegetation, …]. In nature, the overall response of soil water storage cannot be determined by simply observing the response of one factor at a time and subsequently adding the individual observational results together.

Earlier efforts in characterizing the soil spatial variability mainly focused on the spatial similarity in soil properties over a given area using soil classification based on soil survey and conventional statistics [18]. These methods assumed the variation to be random and spatially independent and did not quantify the variability of soil properties with respect to their spatial arrangement, spatial similarity or periodicity [2]. Geostatistical analysis, based on the theory of regionalised variables [19], has been used to characterize spatial similarity in soil properties [3, 20]. The spatial structure or the similarity information as a function of separation distance or scale [11] helps identify autocorrelation in replicating samples, reveals patterns in the data series and identifies the scale of major ongoing processes [21]. However, a necessary assumption in calculating a meaningful variogram, a cornerstone of geostatistical analysis, is that the variable is spatially stationary and the sum of squared differences depends only on the separation of measurements and not on their absolute locations [2, 22].

Different processes controlling the variability in soil properties can contribute differently. Variance contribution from individual processes towards the total variance can be evaluated by transforming the spatial domain information to the frequency domain. The processes that operate at a very fine scale have high spatial frequencies, whereas processes that operate at very broad scales have low spatial frequencies[23]. The contribution of different processes at different scales and their repeated behaviour can be quantified using spectral analysis [12, 24]. Spectral analysis approximates a spatial series by a sum of sine and cosine functions. Each of the functions has an amplitude and a frequency or period. The squared amplitude at a given frequency is equal to the variance contribution of the frequency component to the total variance in the spatial series [25-27]. It deals with the global information or the mean state and cannot examine localized variations or long-term trends [23]. Moreover, spectral analysis assumes that a spatial series is second-order stationary (i.e. the mean and variance of the series are finite and constant). This assumption is generally stricter than the intrinsic stationarity assumption of geostatistics.

Wavelet transformation[28] has been used to examine both long term trends and localized features in soil spatial variation. It enables analysis of multi-scale stationary and/or nonstationary soil spatial variation over a finite spatial domain and has become popular for examining scale and location dependent soil spatial variation [23, 29-33]. A review of the application of the wavelet transform in soil science can be found in Biswas and Si [34]. This method has been extremely useful in examining nonstationary soil spatial variation. However, in examining nonlinear soil spatial variation, the appropriateness of the wavelet transform has been questioned. The Hilbert-Huang transform was introduced into soil science to examine the nonstationary and nonlinear soil spatial variation together [32, 35]. These methods have been useful in examining soil spatial variability. However, they only deal with how the second moment of a variable changes with scales or frequencies.

For normally distributed variables, the second moment plus the average provide a complete description of the variability in the spatial series. The inherent soil variability and the extrinsic factors can cause orders of magnitude of variation in measured soil property values. The presence of intermittent high and low data values often result in distributions deviated from the normal. For these type of distributions (e.g., left skewed distribution), higher moments are needed for a complete description of the variability in the measured property. If we define the q^th moment of Zas$$⟨z^q⟩$$, then when q is positive, the q^th moments magnify the effect of larger numbers in the spatial series Z and diminish the effect of smaller numbers in the data series. On the other hand, when q is negative, the q^th moments magnify the effect of small numbers and diminish the effect of large numbers in the spatial series Z. In this way, by varying the order of the moments, we can look at the magnitude of data values and visualize a better picture of the data series. There is a need to summarize how these moments change with scales or the scaling properties of these moments to compare and simulate spatially-variable soil properties.

Soil properties sometimes vary in space in an irregular manner or exhibit no deterministic patterns. If the irregularity in a variable’s distribution remains statistically similar at all studied scales, the variable is assumed to be self-similar [36]. Self-similarity is closely associated with the transfer of information from one scale to another (scaling). Exploring self-similarity or inherent differences in the scaling properties is important in order to understand the nature of the spatial variability in soil properties. Fractal theory originated by Mandelbrot [37] can be used to investigate and quantitatively characterize spatial variability over a range of spatial scales. Mandelbrot’s early work in the area of geophysics, specifically the characterization of coast lines, showed that the patterns observed at different scales could be related to each other by a power function, whose exponent was called as the fractal dimension. The fractal geometry offered both descriptive and predictive opportunities in the field of soil science [38]. It has been useful in providing a unique quantitative framework for integrating soil biological, chemical and physical phenomena over a range of spatial and temporal scales.

Most of the fractal theory applications in soil science have used a single fractal dimension to characterize the spatial variability over a range of scales [4, 38]. This is known as a monofractal approach - it assumes that the soil spatial distribution can be uniquely characterized by a single fractal dimension. However, the spatial variability in soil properties is the result of the individual or combined influence of soil physical, chemical and biological processes operating at different intensities and scales [11]. Therefore, monofractal distributions are not likely to be prevalent in the landscape [39]. A single fractal dimension will not always be sufficient to represent complex and heterogeneous behavior of soil spatial variations. An extension of the monofractal approach was introduced in soil scienceto describe data with a set of fractal dimensions instead of a single value. A spectrum can be prepared by combining all the fractal dimensions and is known as the multifractal spectrum. The method of characterizing variability based on the multifractal spectrum is known as multifractal analysis [40]. Multifractal behavior is associated with a system where the underlying physics are governed by a random multiplicative process (i.e., successive division of a measure and its spatial support based on a given rule). Therefore, multifractal behavior implies that a statistically self-similar measure can be represented as a combination of interwoven fractal dimensions with corresponding scaling exponents. The multifractal parameters are generally independent of the size of the studied objects [41] and do not assume any specific distribution in the data [42]. Multifractal analysis can transform irregular data into a compact form and amplify small differences among the variables [43, 44]. It uses a wider range of statistical moments, providing a much deeper insight into the data variability structure compared to the methods that use only first two statistical moments. Multifractal analysis can thereforebe used to characterize the variability and heterogeneity in soil properties over a range of spatial scales [44-47]. While the multifractal analysis characterizes the spatial variability in a variable, joint multifractal analysis characterizes the joint distribution of two variables along a common spatial support. It can provide information on the relationships between two variables across different spatial scales [47, 48]. The objective of this chapter is to demonstrate what multifractal and joint multifractal analysis can do in dealing with spatial data series from the published soil science literature. In the next section we briefly describe multifractal and joint multifractal analysis methodology including the working steps and then proceed to describe some applications of these methods in soil science. Finally we close the chapter with a discussion on future prospects for multifractal and joint multifractal methods in soil sciencewithin a brief Conclusions Section.

2. Methods

2.1. Multifractal analysis

Multifractal analysis can be used to characterize the scaling property of a soil variable measured in a direction (such asa transect) as the mass (or value) distribution of a statistical measure on a geometric or spatial support. The geometric or spatial support indicates the extent of the sampling. This can be achieved by dividing the length of the transect into smaller and smaller segments based on a rule that generates a self similar segmentation. One such rule is the binomial multiplicative cascade [36] that can divide a unit interval associated with a unit mass, M (a normalized probability distribution of a variable or a measured distribution as used in a generalized case) into two segments of equal length. This also partitions the associated mass into two fractions, h×M and (1-h)×M and assign them as the left and right segments, respectively. The parameter h is a random variable values between 0 and 1. The new subsets and its associated mass are successively divided into two parts following the same rule. The differences among the subsets are identified using a wide range of statistical moments, which can then be used to determine the multifractal spectra of the measures [36, 48]. Therefore, the multifractal analysis focuses on how the measure of mass varies with box size or scale and provides physical insights at all scales without any ad hoc parameterization or homogeneity assumptions [49]. A detailed description of the multifractal theory can be found in reference[36, 50, and 51]. In this section, for brevity, we will only summarize the computational techniques and concepts commonly used in examining the soil spatial variation.

Past research has indicated that certain properties of a spatial series decrease with increase in scale, following a power law relationship. For example, when all or part of the variogram follows a power law equation of the form $$\gamma \left(h\right)\approx h^{-\zeta }$$the data are scaling in that range—i.e., there is certain degree of statistical scale-invariance [52].

However, the semivariogram only measures the scaling properties of the second moment. Similarly, we can do the same thing for the higher moments such as third, fourth, and so on. Will the scaling properties be the same at higher moments or change with the order of the moments, q? If the scaling properties do not change with q, then we say the spatial series is monofractal, i.e., it only requires a single scaling coefficient to transfer information from one scale to another. If the scaling coefficient changes with q, then the spatial series is multifractal, i.e., it requires multiple scaling coefficients for transferring information over scales. For a spatial series, the scale-invariant mass exponent (structure function),$$\tau (q)$$ , is defined as[46]:

$$⟨{\left[\Delta Z(x)\right]}^q⟩\propto x^{\tau (q)}$$E1

where the symbol “$$\propto$$” indicates proportionality, Z is the spatial series, and x is the lag distance. If the plot of $\tau \left(q\right)$vs.$$q$$ [or $\tau \left(q\right)$curve; Fig. 1a] has a single slope (i.e., a straight line), then the series is a simple scaling (monofractal) type. If the$\tau \left(q\right)$ curve is nonlinear and convex (facing downward), then the series is a multiscaling (multifractal) type (Fig. 1b).

The type of scaling can be examined from the degree of fractality (i.e., monofractal or multifractal) by comparing the$\tau \left(q\right)$curve with a reference curve or a theoretical model (Fig. 1a). One such reference curve (similar to the monofractal type of scaling) or the theoretical model was proposed by Schertzer and Lovejoy [49], which is known as the universal multifractal model or UM model. The UM model simulates an empirical moment scaling function of a cascade process assuming the conservation of mean value. The UM model can be used to compare and characterize the observed scaling properties as reference to the monofractal behavior or scaling. The similarity/dissimilarity can be examined from the goodness-of-fit between the $$\tau (q)$$curve and the UM model using the chi-square test. Statistically significant difference between the curves indicates multifractal scaling and non-significant difference indicates monofractal scaling. The degree of monofractal/multifractal can also be examined by calculating the deviation of the$$\tau (q)$$ curve from the UM model [47, 53]. Large sum of squared difference between the curves indicates multifractal behavior and small sum of squared difference indicates monofractal behavior. The slopes of the regression line fitted to the$$\tau (q)$$ curve (also referred as single fit) can be compared to the slopes of the UM model. Significant difference in the slope indicates multifractal behavior. The$$\tau (q)$$can even be fitted to two regression lines; one for q>0 and another q<0 (also referred as segmented fit). The difference between the means of slopes from segmented fits (for positive and negative q values) can be tested using the Student’s t test. Significant difference between the slopes indicates nonlinearity in the$$\tau (q)$$curve and thus multifractal behavior [47, 53].

Figure 1.

A typical of a) mass exponent graph (or $$\tau (q)$$curve), b) multifractal spectrum (or$$f\left(q\right)$$curve) and c) generalized dimension graph (or Dq curve) used to identify the monofractal and multifractal behavior in any property.

Research has indicated that the q^th order normalized probability measures of a variable (also known as the partition function),$$\mu (q,\epsilon )$$ , vary with the scale size, $$\epsilon$$in a manner similar to Eq. (1)[36, 51], i.e.,

$${\mu }_i(q,\epsilon )=\frac{{\left[p_i(\epsilon )\right]}^q}{{\displaystyle \sum _i{\left[p_i(\epsilon )\right]}^q}}\propto {(\epsilon /L)}^{\tau (q)}$$E2

where$$p_i(\epsilon )$$is the probability of a measure in the $$i^{th}$$segment of size $$\epsilon$$ units, calculated by dividing the value of the variable in the segment by the whole support length of L units. In simpler terms, $$p_i(\epsilon )$$measures the concentration of a variable of interest (e.g., clay content, organic carbon, etc.) in a given segment relative to the whole support length. The$$\tau (q)$$function in Eq. (2) is given a new name, ‘the mass exponent,’ because it relates the probability of mass distribution in a given segment to the size of the segment (scale) and is used widely in multifractal analysis.

The multifractal spectrum, $$f\left(q\right)$$(Fig. 1b), which is the fractal dimension of the subsets of segments of size $\epsilon$units with a coarse Hölder exponent (local scaling indices) of $\alpha$if calculated to the limit as$$\epsilon \to 0$$, can be expressed as [36, 50],

$$f(q)=\underset{\epsilon \to 0}{{\displaystyle \lim }}{\left(\log \left(\frac{\epsilon }{L}\right)\right)}^{-1}{\displaystyle \sum _i{\mu }_i(q,\epsilon )\log {\mu }_i(q,\epsilon )}$$E3

and the local scaling indices, $\alpha$, are given by

$$\alpha (q)=\underset{\epsilon \to 0}{{\displaystyle \lim }}{\left(\log \left(\frac{\epsilon }{L}\right)\right)}^{-1}{\displaystyle \sum _i{\mu }_i(q,\epsilon )\log p_i(\epsilon )}$$E4

Note that f(α) can also be determined through the Legendre transform of the τ(q) curve:

$$f(\alpha (q))=q\alpha (q)-\tau (q)$$E5

The multifractal spectrum is a powerful tool in portraying the variability in the scaling properties of the measures (e.g., clay content, organic carbon, etc.). The spectrum also enables us to examine the local scaling property of soil variable. The width of the spectrum (α_max - α_min) is used to examine the heterogeneity in the local scaling indices (Fig. 1b). The wider the spectrum, higher the heterogeneity in the distribution of the soil variable. Similarly, the height of the spectrum corresponds to the dimension of the scaling indices. Small f(q) values correspond to rare events (extreme values in the distribution), whereas the largest value is the capacity dimension that is obtained at q = 0.

For many practical applications indicator parameters are selected and used, in addition to the multifractal spectrum (f(q) vs. $\alpha$(q)), to describe the scaling property and variability of a process. Generalized dimensions,$$D_q$$ , (Fig. 1c) are often used to provide indicator parameters. The$$D_q$$of a multifractal measure is calculated as

$$D_q=\frac{1}{q-1}\underset{\epsilon \to 0}{{\displaystyle \lim }}\frac{\log {\displaystyle \sum _i{p}_i}(\epsilon )}{\log (\epsilon )}$$E6

The D_q value at q = 0, D₀, is called the capacity dimension or the box counting dimension of the geometric support of the measure. The value atq = 1,$$D_1$$ , is referred to as the information dimension and provides information about the degree of heterogeneity in the distribution of the measure – this is analogous to the entropy of an open system in thermodynamics [54]. Sometimes,$$D_1$$ is also known as the entropy dimension. A value of D₁ close to unity indicates the evenness of measures over the sets of a given cell size, whereas a value close to 0 indicates a subset of scale in which the irregularities are concentrated. The$$D_2$$, known as the correlation dimension, is associated with the correlation function and measures the average distribution density of the measure [55]. For a distribution, with simple scaling (monofractal), the D₁ and D₂ become similar to the D₀, the capacity dimension. The value of D₀ = D₁ = D₂ indicates that the distribution exhibit perfect self-similarity and is homogeneous in nature (Fig. 1c). However, for multifractal type scaling the order generally becomes D₀>D₁>D₂. The D₁/D₀ value can also be used to describe the heterogeneity in the distribution [56]. A D₁/D₀ value equal to 1 indicates exact monoscaling of the distribution, which means that all fractions take equal values at different scales.

2.1.1. Multifractal analysissteps

1. Calculate the probability measure (p) from a linear distance for a transect, or from a rectangle for an area. Regarding the minimum number of samples required to carry on the analysis, the multifractal analysis method has the flexibility over other methods such as Fourier transform or wavelet transform, which generally require a larger dataset with regular interval between samples.

2. For each n = 2, 4, 8, 16,... until the unit is not dividable. To calculate the probability measure (p), sum up the values of all the points in the segment, and divide by the total of all segments $$\left(p_i(\epsilon )=\frac{m(\epsilon ,i)}{M}\right)$$ where i is the index for i^th segment when the whole transect is divided into n segments with a unit length =$$\epsilon \left(i=1,\cdots ,n\right)$$; m(ε, i) is the sum of measurements at all points in the i^th segment of length ε; M is the total of all points along the transect.

3. For certain q values or statistical moments (e.g., -20 to 20, which are selected based on the nature of the data to be analyzed), calculate μ_i(q, ε) using Eq. (2).

4. Calculate $$Ff(q,\epsilon )={\displaystyle \sum _i{\mu }_i(q,\epsilon )\log {\mu }_i(q,\epsilon )}$$

5. Calculate $$F\alpha (q,\epsilon )={\displaystyle \sum _i{\mu }_i(q,\epsilon )\log p_i(\epsilon )}$$

6. Check if the partition function μ(q, ε) obeys the power law or if the log-log plot of the partition function and distance is linear. If they are not, then they are not multifractal and no further multifractal analysis is needed. If they are, continue the following.

7. Calculate τ(q) as the intercept of linear regression of $$\log \left[{\mu }_i(q,\epsilon )\right]$$vs.$$\log (\epsilon ).$$

8. Calculate f(q) as the intercept of linear regression of $$Ff(q,\epsilon )$$ vs. $$\log (\epsilon ).$$

9. Calculate α(q) as the intercept of linear regression of $$F\alpha (q,\epsilon )$$ vs. $$\log (\epsilon ).$$

10. Calculate D_q as the intercept of the linear regression of $$FD(q,\epsilon )$$ vs. $$\log (\epsilon ).$$

11. Plot τ(q) as a function q, f(q) as a function of α(q), and D_q as a function of q.

All of these calculations can be implemented in MathCad, MATLAB, or SAS.

2.2. Joint multifractal analysis

While the multifractal analysis characterizes the distribution of a single variable along its spatial support, the joint multifractal analysis can be used to characterize the joint distribution of two or more variables along a common spatial support. Similar to the multifractal analysis, the length of the datasets (for example, 2 datasets) is divided into a number of segments of size ε. The probability of the measure of the i^th segment of the first variable is$$P_i\left(\epsilon \right)\infty {\left(\epsilon /L\right)}^{\alpha }$$and for the second variableis$$R_i\left(\epsilon \right)\infty {\left(\epsilon /L\right)}^{\beta }$$, where α and β are the local singularity strength corresponding to $$P_i\left(\epsilon \right)$$ and$$R_i\left(\epsilon \right)$$, respectively. The partition function (the normalized μ measures) for the joint distribution of $$P_i\left(\epsilon \right)$$and$$R_i\left(\epsilon \right)$$, weighted by the real numbers q and t can be calculated as [47, 50, 51];

$${\mu }_i(q,t,\epsilon )=\frac{p_i{(\epsilon )}^q\cdot r_i{(\epsilon )}^t}{{\displaystyle \sum _{j=1}^{N(\epsilon )}\left[p_j{(\epsilon )}^q\cdot r_j{(\epsilon )}^t\right]}}$$E7

The local scaling indices (coarse Hölder exponents) with respect to two probability measures$$P_i(\epsilon )$$ and$$R_i(\epsilon )$$, which are representedby$$\alpha (q,t)$$and$$\beta (q,t)$$ respectively, are calculated as follows:

$$\alpha (q,t)=-{\left[\log (N(\epsilon ))\right]}^{-1}{\displaystyle \sum _{i=1}^{N(\epsilon )}\left[{\mu }_i(q,t,\epsilon )\cdot \log (p_i(\epsilon ))\right]}$$E8

$$\beta (q,t)=-{\left[\log (N(\epsilon ))\right]}^{-1}{\displaystyle \sum _{i=1}^{N(\epsilon )}\left[{\mu }_i(q,t,\epsilon )\cdot \log (r_i(\epsilon ))\right]}$$E9

The dimension (i.e.,$$f(\alpha ,\beta )$$) of the set on which $$\alpha (q,t)$$ and $$\beta (q,t)$$are the mean local exponents of both measures is given by

$$f(\alpha ,\beta )=-{\left[\log (N(\epsilon ))\right]}^{-1}{\displaystyle \sum _{i=1}^{N(\epsilon )}\left[{\mu }_i(q,t,\epsilon )\cdot \log ({\mu }_i(q,t,\epsilon ))\right]}$$E10

When q or t is set to zero, the joint partition function shown in Eq. [7] reduces to the partition function of a single variable, and hence the joint multifractal spectrum defined by Eq. [10] becomes the spectrum of a single variable. When both q and t are set to zero, the maximum $f\left(\alpha ,\beta \right)$is attained, which is the dimension if all the segments contain the same concentration of mass. Different pairs of $\alpha$and $\beta$can be scanned by varying the parameters q and t. Because, high q or t values magnify large values in the data and negative q or t values magnify small values in the data, by varying q or t, we can examine the distribution of high or low values (different intensity levels) of one variable with respect to that of the other variable. Pearson correlation analysis can be used to quantitatively illustrate the variation of the scaling exponents of one variable with respect to another variable across similar moment orders. Because $f\left(\alpha , \beta \right)$represents the frequency of the occurrence of a certain value of $\alpha$and a certain value of$\beta$, high values of $f\left(\alpha , \beta \right)$signifies a strong association between the values of $\alpha$and the values of$\beta$. By perturbingq and t, we can examine the association of similar values (highs vs. highs or lows vs. lows) of $\alpha$and $\beta$as well as dissimilar values (highs vs. lows) of $\alpha$and$\beta$. Generally, a contour graph is used to represent the joint dimensions, $f\left(\alpha , \beta \right)$of the pair of variables (Fig. 2). The bottom left part of the contour graph shows the joint dimension of the high data values of the two variables, while the top right part represents the low data values [57]. The diagonal contour with low stretch indicates strong correlation between values corresponding to the variables in the vertical and horizontal axis (Fig. 2).

Figure 2.

A typical of multifractal spectra for joint distribution of two variables with a) strong correlation and b) weak correlation.

2.2.1. Joint multifractal analysissteps

1. Calculate the probability measure (p) for variable Y from a linear distance for a transect or a rectangle from an area.

2. For each n = 2, 4, 8, 16,... until the unit is not dividable, calculate the probability measure (p), sum up values of all the points in the segment, and divide by the total of all segments $$\left(p_i(\epsilon )=\frac{m_Y(\epsilon ,i)}{M_Y}\right)$$ where i is the index for i^th segment when the whole transect is divided into n segments with a unit length =$$\epsilon \left(i=1,\cdots ,n\right)$$; Similarly, we can calculate the probability measure for the variable Z:$$r_i(\epsilon )=\frac{m_Z(i)}{M_Z}$$.

3. For certain q values (e.g., -20 to 20) and t values (e.g., -20 to 20), calculate μ_i(q, t, ε) using Eq. (7).

4. Calculate $$Ff(q,t,\epsilon )={\displaystyle \sum _i{\mu }_i(q,t,\epsilon )\log {\mu }_i(q,t,\epsilon )}$$

5. Calculate $$F\alpha (q,t,\epsilon )={\displaystyle \sum _i{\mu }_i(q,t,\epsilon )\log p_i(\epsilon )}$$

6. Calculate $$F\beta (q,t,\epsilon )={\displaystyle \sum _i{\mu }_i(q,t,\epsilon )\log r_i(\epsilon )}$$

7. Calculate τ(q, t) as the intercept of linear regression of $$\log \left[\mu (q,t,\epsilon )\right]$$vs.$$\log (\epsilon )$$ (Eq. (7)).

8. Calculate f(q, t) as the intercept of linear regression of $$Ff(q,t,\epsilon )$$ vs. $$\log (\epsilon )$$(Eq. (8)).

9. Calculate α(q, t) as the intercept of linear regression of $$F\alpha (q,t,\epsilon )$$ vs. $$\log (\epsilon )$$(Eq. (9)).

10. Calculate β(q, t) as the intercept of linear regression of $$F\beta (q,t,\epsilon )$$ vs. $$\log (\epsilon )$$(Eq. (10)).

11. Plot τ(q) as a function q and contour plot of f(α, β).

Again, all these calculations can be implemented in MathCad, MATLAB, or SAS.

2.3. Comments on multifractal and joint multifractal analysis

• Multifractal analysis for two-dimensional or three dimensional fields is the same as for a transect. However, the support ε becomes an area or volume.

• Multifractal analysis, like spectral analysis, is based on the global statistical properties of spatial series. Therefore, localized information is lost, which is different from wavelet analysis. However, multifractal analysis provides information regarding the higher moments and how the higher moments change with scale.

• Multifractal analysis does not require regularly spaced samples. Any sampling scheme can be analyzed by multifractal analysis.

• For highly spatially variable fields, the probability p for some locations may be very small or even zero. As such, the negative power of p can be very large and the partition function will be dominated by this single value. In this case, the multifractal method may diverge and the process of division needs to be stopped.

• Joint multifractal analysis can be used for the simultaneous analysis of several multifractal measures existing on the same geometric support, and hence for quantifying the relationships between the measurements studied. Joint multifractal analysis is based on the assumption that the individual variable is multifractal.

• Simulation of a synthetic field according to the measured multifractal distribution may enhance our understanding of the effects, spatial scaling, and spatial variability of soil properties on various soil processes.

3. Application of multifractal and joint multifractal analysis in soil science

Fractal theory [37] has been used to investigate and quantitatively characterize spatial variability over a large range of measurement scales in different fields of geosciences including soil science [4, 58]. A detailed review of the applications of fractal theory in soil science can be found in reference[38]. The fractal theory applications in soil science used monofractal approach, which assumes that the soil spatial distribution can be uniquely characterized by a single fractal dimension. However, a single fractal dimension might not be sufficient to represent complex and heterogeneous behavior of soil spatial variations. Multifractal analysis, which uses a set of fractal dimensions instead of one, is useful in characterizing complex soil spatial variability. Use of multifractal analysis in explaining soil spatial variability started during the late 1990’s but only become popular in soil science after 2000. Among the earlier works, Folorunso et al. [59] used multifractal theory for analyzing spatial distribution of soil surface strength and Muller [60] used multifractal analysis for characterizing pore space in chalk materials. Since then, multifractal analysis has been used in soil science to study various issues including spatial variability of soil properties (e.g., physical, chemical, hydraulic), soil groups and pedotaxa, soil particle size distribution, soil surface roughness and microtopography, effect of tillage activities, crop yields, soil porosity and pore size distribution, flow and transport of water and chemical in soil, infiltration, and downscaling soil water information from satellite images. While multifractal analysis is used to characterize the spatial variability of soil properties over a range of scales, the joint multifractal analysis can be used in soil science to characterize the joint distributions between soil properties over a range of scales. In this section we review previous applications of multifractal and joint multifractal analysis in soil science.

3.1. Multifractal analysis

Various soil properties have shown multifractal behavior and the variability in those soil properties has been characterized using multifractal analysis. For example, Kravchenko et al. [44] was one of the first to report the multifractal nature of soil-test phosphorus (P), potassium (K), exchangeable calcium (Ca), magnesium (Mg) and cation exchange capacity (CEC). The study used a 259 ha agricultural field in central Illinois, USA. One set of grid size was used and the authors reported only one unique set of multifractal spectra for each soil property. Multifractal parameters were studied over a range of moment orders (q) from -15 to 15. The minimum value of multifractal parameters $\alpha$(scaling indices) and $f\left(\alpha \right)$(multifractal spectra), i.e., ${\alpha }_{min}$and $f\left({\alpha }_{min}\right),$corresponded to q = 15, and the maximum values, ${\alpha }_{max}$and $f\left({\alpha }_{max}\right)$corresponded to q = -15. A high value of ${\alpha }_{max}-{\alpha }_{min}$ indicated a wider opening of the multifractal spectrum and thus the multifractal nature of the soil properties except organic matter (OM) content and soil pH (Fig. 3) [44]. A very small value of ${\alpha }_{max}-{\alpha }_{min}$ indicated the monofractal nature of OM and pH (Fig. 3). The generalized fractal dimensions were also calculated for all positive q values. An excellent fit between the theoretical fractal dimension model and the D(q) curve also indicated the multifractal nature of the soil properties except OM and soil pH. A high correlation between the multifractal parameters and exploratory statistics (e.g., coefficient of variations, skewness, kurtosis) or geostatistical parameters (e.g., nugget, range) provided comprehensive information on the major aspects of data variability [44]. Multifractal spectra of the soil properties studied carried a large amount of spatial information and allowed quantitative differentiation between the soil variability patterns. Kravchenko et al. [44] highlighted the usefulness of multifractal parameters of soil properties in the interpolation and mapping of those properties. Interpolation methods based on the multifractal scaling relationship improves the mapping of soil properties.

Figure 3.

Multifractal spectra, f(q), for a) phosphorus, b)potassium, c) organic matter, d) pH, e) exchangeable Ca, f) exchangeable Mg, and g) cation exchange capacity (adopted from 44).

Various other studies from different parts of the world also indicate the multifractal behavior of soil properties. Caniego et al. [61] studied spatial variability in soil electrical conductivity (EC), pH, OM content and depth of Bt horizon along three transects - two short transects (each 33 m long) with high intensity sampling (sampling interval 25 cm), and a long transect (3 km) with 40 m sampling interval. Though the authors reported the multifractal nature of the soil properties along the long transect, the variability along the short transects were more homogeneous. The variability along the short transect was explained by simple random fractal noise (close to a white noise) and thus a single scaling index was consideredsufficient to transform information from one scale to another [61]. The variability in the long transect might have a deterministic component reflecting changing geological features and differences in soil forming factors with distance. The presence of nonlinearity together with the environmental features along the long transect resulted in the observed multifractal behavior of soil properties. Monofractal behavior of sand and silt was reported by Wang et al. [62]. Zeleke and Si [52, 57] characterized the distribution of bulk density (Db), sand content (SA) and to some extent silt (SI) content as monofractal along a gently sloping 384 m long transect in semi-arid central Saskatchewan, Canada. However, these authors did report multifractal behavior of clay (CL), saturated hydraulic conductivity (K_s) and organic carbon (OC) from the same study. For example, Fig. 4 shows the $\tau \left(q\right)$(mass exponents) and $f\left(q\right)$(multifractal spectra) curve of the soil variables studied by Zeleke and Si [52]. The linearity of the τ(q) curve of Db, SA and SI indicated a monofractal nature, while the convex shape of the τ(q) curve of CL, OC, and K_s indicated multifractal behavior. Similarly, a large value for${\alpha }_{max}-{\alpha }_{min}$ of CL, OC, and K_s, and thus the wide opening of the curve, indicated the multifractal nature of these properties (Fig. 4) [52]. Wang et al. [62] also reported highly multifractal behavior of OC and CL content along a transect and monofractal behavior of SA and SI. A decrease in the variability of EC was reported with an increase in the EC value itself [63].

Figure 4.

a) Mass exponents of different variables. The solid line is the UM model passing through$\tau \left(0\right)$. b) The multifractal spectra of different variables. The q = -6 to 6 at 0.5 increments at both cases. Db is the bulk density, SA is the sand content, SI is the silt content, CL is the clay content, OC is the organic carbon content, Ks is the saturated hydraulic conductivity. (adopted from 52)

The multifractal behavior of effective saturated hydraulic conductivity in a two-dimensional spatial field was reported by Koirala et al. [64] from a numerical simulation study. The authors used a normalized mass fraction of a randomized multifractal Sierpinski carpet to represent the areal distribution of saturated hydraulic conductivities in an aquifer. The effective saturated hydraulic conductivity was related to the generalized dimensions of the multifractal field. This relationship may be helpful to predict the effective saturated hydraulic conductivity of an aquifer from an empirical Dq spectra determined by the multifractal analysis [64]. The multifractal nature of the hydraulic conductivity was also reported by Liu and Molz [46].

The multifractal behavior of soil water content [53, 63,65, 66], soil water retention at different suction [57,62, 67,68] and theoretical water retention parameters (e.g., van Genuchten$\alpha$ and n) [62, 67-71] was also reported in literature. The volumetric as well as gravimetric soil water content was found to exhibit multifractal behavior [65]. The variability in soil water content decreases with the increase in soil water and increases with the increase in the representative area of the measurement [66]. Though the water retention at saturation showed monofractal behavior, water retention at -30 kPa and -1500 kPa showed multifractal behavior [57, 67]. Filling up the pore with water during saturation irrespective of the size of the pores might result in characteristic that is persistent over a range of scales. However, at -30 kPa and -1500 kPa, the relation between pore size and pore water retention might result in multifractal behavior [57, 67]. While the van Genuchten [72] water retention parameter n (considered as relating to pore size distribution) showed monofractal behavior, the parameter $\alpha$(considered as the inverse of air entry potential) showed weak multifractal behavior [62, 69].Similar behavior of water retention parameters was reported by Liu et al. [70, 71]. However, Guo et al. [68] reported multifractal behavior of both water retention parameters,nand$\alpha$.

The fractal behavior of soil water estimates derived from the satellite images were reported earlier [73, 74]. The relationship between the variance in soil water and the area of measurement or aggregation area (or scale) indicated self-similarity or scale invariance in soil water estimates. Later this relationship was used to develop models for downscaling information on soil water estimates from satellite images. A number of studies have developed and used multifractal models to downscale soil water estimates from satellite images to represent small areas[75-79]. Often the passive remote sensing images provide estimates of soil water for a large area (coarse spatial resolution;e.g. 25-50 km). Downscaling models are necessary tools to characterize and reproduce soil water heterogeneity from the remotely sensed estimates. The presence of multifractal behavior helped developing downscaling models [76].

Prediction of soil water storage as affected by soil microtopography or microrelief can also be made using the multifractal approach [80-81]. Soil surface roughness and microtopography showed multifractal behavior, the quantification of which help characterizing spatial features in topography and water storage in micro-depressions. The soil surface roughness created by tillage operations can determine soil strength, penetration of roots and susceptibility to wind and water erosion. The multifractal behavior of soil surface roughness can be used to explain the structural complexity created by tillage operations and the effect of various tillage implements. Various studies have reported the multifractal behavior of soil surface roughness [80-85]. Generally, these studies measured the soil surface roughness in a very small area. The degree of multifractality was found to increase at higher statistical moments [84]. The variability in the distribution of soil units also showed multifractal behavior [86-87]. Information on the variability within and between soil pedological units (pedotaxa) shows promise in analyzing and characterizing the complexity of soil development at multiple scales.

A large number of studies reported the multifractal behavior of soil particle size distribution [56, 88- 95]. Most of the studies used soil particle size distributions measured using laser diffraction. One of the studies also used a piece-wise fractal model for very fine and coarse particle sizes [95]. Multifractal behavior of soil particles sizes at different land can be used to characterize soil physical health and its quality [90].

Distribution of soil particles size determines the porosity, which in turn determines the flow and transport of water and chemicals in soil. Soil thin sections have been used to study the pore arrangements and distributions in soil [96-98]. The pore geometry showed multifractal behavior, which is useful for classification of soil structure and determining the fluid flow parameters [96]. The multifractal nature of soil porosity was used to identify the representative elementary area using photographs of soil and confocal microscopy [99]. The use of two dimensional binary or grayscale images[96-98] of soil thin sections helped in characterizing the pore structure and movement of water in soil. With the advancement of technology, the images from Magnetic Resonance Imaging (MRI) [102] or Computer Tomography (CT) [101, 102] helped in characterizing the soil as a porous media and the fluid flow through it. Three-dimensional images of soil systems helped to identify soil pores and their connectivity in three-dimensions, which in turn helped understand the movement of water in soil and the characterization of preferential flow paths [100]. The multifractal behavior of soil pore systems is often used to develop models for creating simulated media, which helps to study flow and transport of water and chemicals in soil or other porous media [103].

The infiltration of water into soil often shows multifractal behavior. Use of dye is a common approach to study the preferential flow paths of water in soil. The distribution of dye along with the infiltration water helped illustrate multifractal behaviour[104-106]. The multifractal nature of infiltration helped in identifying the input parameters for various rainfall-runoff models [107-109]. Perfect et al. [109] used a multifractal model to develop a simulated media for upscaling effective saturated hydraulic conductivity.

The radioactive cesium (¹³⁷Cs) fallout at small spatial scales has been found to be multifractal in nature [111]. Grubish [111] reviewed research on the ¹³⁷Cs fallout and reported that the multifractal nature stemmed from the erosion and deposition of soil materials, which showed a log-normal distribution. The behavior of soil chemical processes such as nitrogen absorption isotherms was also found to be multifractal [110]. The authors studied the isotherm characteristics from 19 soil profiles in a tropical climate. The asymmetric singularity spectra (shifted to the left) indicated the highly heterogeneous and anisotropic distribution of the measure (Fig. 5). The generalized dimension spectra (Dq) was used to discriminate among soil groups with contrasting properties, arising from different pedological processes such as Ferralsols (Latosols) and soil with argillic (Bt) horizons and/or an abrupt textural change.

Various studies also reported the multifractal behavior of soil landscape properties. For example, Wang et al. [112] found the multifractal nature of the spatial and temporal patterns of land uses in China. The multifractal pattern was also found in different terrain indices. While the upslope area and wetness index were multifractal, the relative elevation was monofractal [47, 113]. This may be due to inclusion of multi-scale characteristics in calculating secondary terrain indices such as the wetness index. The spatial variability of crop yield also found to exhibit multifractal in nature [47, 113-115]. A stochastic simulation study showed that the spatial variability in the production of wheat crop was multifractal, while the production of corn was monofractal [115].

Figure 5.

A multifractal spectrum generated using duplicate samples from two soil horizons that exhibit the lowest and the highest mean Hölder exponents of order zero (adopted from 110).

Figure 6.

Multifractal spectra of the joint distribution of saturated hydraulic conductivity (vertical axis) and other soil properties (horizontal). Contour lines show the joint scaling dimensions of the variables. Db is the bulk density, SA is the sand content, SI is the silt content, CL is the clay content, OC is the organic carbon content, Ks is the saturated hydraulic conductivity (adopted from 52)

3.2. Application of joint multifractal analysis in soil science

Joint multifractal analysis has been used to characterize the joint distribution between two soil properties over a range of scales. Often soil hydraulic properties are predicted from more easily measured soil physical properties using pedotransfer functions. The scale dependence of soil physical properties often creates issues in the prediction. A large body of literature used joint multifractal analysis to characterize the joint distribution between soil physical and hydraulic properties [52, 57, 62, 63, 66, 68-71]. Zeleke and Si [52] studied the relationship between the saturated hydraulic conductivity and various soil physical properties. Figure 5 shows the contour lines for the joint scaling dimensions of K_s and Db. There appear to be some relationships between the scaling dimensions of K_s and Db for both high and low data values, which is evident from the slightly diagonal feature of the plots and the high correlation coefficients between the scaling indices of the two variables (Fig. 6). The highest correlation coefficient (r = -0.57) between scaling indices of Db and K_s was obtained for the high data values of the two variables. However, for joint multifractal spectra of SA and K_s, the contour line is nearly concentric, indicating there are no preferred associations between values of exponents for K_s and sand content. Therefore in this case, prediction of the scaling properties of K_s from those of sand content would not be recommended [52].

Figure 7.

Multifractal spectra of the joint distribution of van Genuchten $\alpha$(vertical axis) and the soil texture (horizental axis). (adopted from 69).

Various theoretical hydraulic parameters (e.g., van Genuchten $\alpha$and n) are used to characterize soil water retention capacity. Often these parameters are predicted from soil physical properties potentially introducing issues of scale dependency. The joint distribution between theoretical soil hydraulic parameters and various soil physical properties were studied using joint multifractal analysis [69]. The contour plot for van Genuchten $\alpha$and sand content are elliptical and tilted to the left(Fig. 7). However, the small difference between the long axis and short axis of the ellipse indicated that the van Genuchten α was not strongly correlated to the sand content over multiple scales. Similarly, there was weak correlation between the van Genuchten α and the silt. However, the elliptical contour tilted to the right indicated a strong positive correlation between the van Genuchten α and clay content (Fig. 7). The small diameter of the contour across the tilt compared with that along the tilt indicated low variability in the joint scaling exponents between the van Genuchten α and clay content. This means that the van Genuchten α and clay content had strong relationship at all intensity levels (high vs. high and low vs. low data values). This may be because the dominant factor of aggregation in loam soil is clay content. Characterizing the relationship at multiple scales using joint multifractal analysis is different from the traditional correlation analysis, which only explains the variability at the scale of measurement. However, variability analyses at a single scale may not accurately reflect the spatial patterns and processes of the studied variables.

The relationship between crop yield and various terrain indices over a range of scales were also studied using joint multifractal analysis [47, 48]. Quantification of the spatial variability in the crop yield and the yield affecting factors has important implications in precision agriculture. Topography is often considered as one of the most important factors affecting yields, and topographical data are much easier to obtain than time and labor-consuming measurements of soil properties. Kravchenko et al. [48] used the joint multifractal analysis to differentiate between yield-distributions corresponding to field locations with high and low slopes. The authors observed larger yields at low slope locations and a wide range of yields at sites with moderate and high slopes during four growing seasons with moderate and dry weather conditions. During the wet growing season, lower yields prevailed at locations with low slopes [48]. The applicability of joint multifractal analysis to study crop yield and topography relationships and characterize spatially distributed data has also been reported by Zeleke and Si [47]. The authors reported that the distribution of yield and biomass was better reflected in upslope length than other topographic indices such as wetness index and relative elevation. This implies that the upslope length can be used as a guideline for varying production inputs such as fertilizer, especially when detailed recommendations at a given scale of interest are not available. The strong relationship between the upslope length and the crop yield implies that one can make any site specific prediction of the final yield (before the harvest) using the upslope length regardless of scale [47]. The variability in crop yield with respect to the nitrate nitrogen level in soil was also studied using joint multifractal analysis [114]. Various authors also developed models to better predict the crop yield based on the multifractal nature of the crop yield and terrain indices [115].

Leaf area index (LAI) is often used as the growth indicator for plants. Banerjee et al. [113] studied the relationship between LAI of pasture and various terrain indices including relative elevation, wetness index, upslope length and distance. There was high correlation between LAI and the wetness index and upslope length over a range of scales. The study indicated that the spatial heterogeneity in LAI is well reflected in the variability of wetness index and upslope length [113] enabling the prediction of variability in LAI and thus the grass production from the variability in terrain indices.

4. Conclusion and future prospects

In this chapter, we reviewed various applications of multifractal and joint multifractal analysis to spatial studies in soil science. Spatial variation of soil has spatial dependence, periodicity, nonstationarity, nonlinearity and many other characteristics. Geostatistics quantifies spatial dependence and spectral analysis can analyze the scale information, but loses spatial information. The wavelet transform can be used to examine the spatial variability of nonstationary series, and the Hilbert-Huang transform can examine the spatial variability of nonstationary and nonlinear series concurrently. These methods can only deal with the second moment of a variable change with scales or frequencies. For a normally distributed variable, the second moment and the average provide a complete description of the variability in a spatial series. However a second order moment can only provide a poor characterization of the variability that occurs in non-normal distributions. For distributions other than normal, we need higher moments for complete characterization of soil spatial variability.

Fractal theory can be used to investigate and quantitatively characterize spatial variability over a large range of measurement scales. According to the fractal theory, the properties that are observed at different scales are related to each other (self-similar) by a power function, whose exponent is called the fractal dimension. When a single fractal dimension is used to characterize soil spatial variability, it is considered as monofractal. However, soil variability can occur in different intensities at different scales and a single fractal dimension may not be sufficient. The multifractal theory implies that a statistically self-similar measure can be represented as a combination of interwoven fractal dimensions with corresponding scaling exponents. A multifractal spectrum is prepared by combining all the fractal dimensions, which can be used to characterize the variability and heterogeneity in a variable over a range of scales. A large number of fractal sets provides in-depth characterization of soil spatial variability compared to that using a single fractal dimension. The multifractal approach is independent of the size of the studied variable and does not require any assumption about the data following any specific distribution. Therefore, the multifractal approach has been used to characterize spatial variability of soil properties, soil groups, soil water, crop yield, and terrain properties. The multifractal nature of the variability has been used to downscale soil water estimates from satellite images. The observed self-similar nature of pore sizes and their distribution has been used to characterize the flow and transport of water and chemicals through soil.

While multifractal theory is used to characterize the variability in any one soil property, the joint multifractal theory was used to characterize the joint distribution of two variables along a common spatial support. Various soil properties are difficult to measure as well as time, cost and labor intensive, while other properties are easy to measure. The variability in one property that can be easily measured can be used to predict the variability in other properties that are difficult to measure. The joint distribution of soil hydraulic properties (e.g., saturated hydraulic conductivity, theoretical hydraulic parameters) with soil physical properties (sand, silt, clay, bulk density, organic carbon) has been used to predict the variability in soil hydraulic properties. Similarly, crop yield has been predicted from its relationship with various terrain indices and other chemical properties.

Multifractal and joint multifractal analyses are versatile tools for characterizing soil spatial variability at multiple scales. Use of the multifractal methods with other approaches such as wavelet analysis would provide extensive information on soil spatial variability at different scales and locations. Multifractal analysis deals with global information and the wavelet transform can deal with local information. Joining the methods would help complete characterization of soil spatial variability [67, 116-118]. Often the localized trends and nonstationarities in the spatial series create a challenge in the scaling analysis. Moreover, nonstationarity can introduce a superficial scaling in the soil properties. In this situation, the intrinsic scaling property from variation in a given soil attributeneeds to be separated from the undue influence of larger scale trends. Use of detrended fluctuation analysis together with multifractal analysis provides an opportunity to characterize the intrinsic variability of soil properties at field scale resulting from the interaction of all underlying processes [53]. Similarly, the use of a multifractal detrended moving-average can be an opportunity to identify the intrinsic scales of variability [119], and multifractal detrended moving-average cross-correlation analysis can be used to analyze and compare the variability between two soil attributes[120].

Acknowledgement

The funding for this work was provided through a CSIRO Land and Water post doctoral fellowship. The thoughtful comments from Dr. Gallant are highly appreciated.