For realization of that research data of the seismic catalogue of Tashtagol mine during two years from June 2006 to June 2008 had been used. As a data set we used the space-time coordinates for all dynamical events-responses of the massif occur in that time period inside the mine field and also explosions, which had been developed for massif outworking and values of energy which had been fixed by the seismological station. In our analysis we divided the whole mine field in two parts (figures 6, 7). The events-responses had been taken into account from horizons – 140 m, – 210 m, – 280 m, – 350 m. According to the catalogue, explosions had been provided on the south-east place – on horizons + 70 m, 0 m, – 70 m, on other places – on all listed higher horizons. The whole catalogue had been

**3. Investigation of non-linear dynamics of rock massive, using**

Ep – absorbed energy (joules), Ev-dissipated energy (joules), Vmax – the maximum dimension of the volume, where had been observed dynamical events.

**Figure 4.** The response of the massif by the man-made influence: 07. 10 – 23. 12 2001, block23, horizon( – 280/ – 210) (I) 19. 01 – 20. 03 2003, block24, horizon ( – 280/ – 210)(II) X-axis– number of explosions

**Figure 5.** The response of the massif by the man-made influence: 27. 02 – 09. 07. 2000, block7, horizon ( – 210/ – 240) (I), 21. 12 – 21. 03 2003 – 2004, block6, horizon ( – 210/ – 140)(II). X-axis – number of explosions The legend of Figure 5 is the same as in Figure 4.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 97

**Figure 8.** Distribution of absorbed (1 – 2) and dissipated (3 – 4) energy of the whole mine field during

**Figure 9.** Distribution of absorbed (1 – 2) and dissipated (3 – 4) energy of the whole mine field during

the period II 13. 01. 2007 – 17. 05. 2008. The horizontal axes is time per days.

the period I 03. 06. 2006 – 13. 01. 2007. The horizontal axes is time per days

**Figure 6.** Plan of horizon – 210, southern place.

**Figure 7.** Plan of horizon – 210, northern place.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 97

96 Fractal Analysis and Chaos in Geosciences

**Figure 6.** Plan of horizon – 210, southern place.

**Figure 7.** Plan of horizon – 210, northern place.

**Figure 8.** Distribution of absorbed (1 – 2) and dissipated (3 – 4) energy of the whole mine field during the period I 03. 06. 2006 – 13. 01. 2007. The horizontal axes is time per days

**Figure 9.** Distribution of absorbed (1 – 2) and dissipated (3 – 4) energy of the whole mine field during the period II 13. 01. 2007 – 17. 05. 2008. The horizontal axes is time per days.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 99

The whole interval of research had been divided into four periods: from 03. 06. 2006 – 13. 01. 2007 (period I), from 14. 01. 2007 – 17. 05. 2008 (period II), from 24. 05. 2008 – 26. 07. 2009 (period III), from 28. 06. 2009 – 18. 07. 2010 (period IV). (Figures 8 – 11).

We would like to note such peculiarities for the fourth period (Figure 11): during the period 100 days, beginning from 100 days from the beginning of the analyzed period and finishing by 200 days, the explosions had been provided in the northern and in the southern part of the mine field approximately of equal intensity, however the energy of massif response in the southern part is significantly larger, than from the northern part. During next 50 days the explosions had been provided in the southern part, but the energy of the massif response in the northern part and in the southern part are approximately equal. During the period from 300 days to 400 days the explosions had been provided mainly in the northern part of the mine field, the distribution of the massif response energy in the northern part practically corresponds to the distribution of the

**Period I**

Ep(joule) Ev(joule) Ep(joule) Ev(joule) R(Ep, Ev)(A) R(Ep, Ev)(B) R(Ev, Ev)

3. 1·107 3. 1·106 2. 04·109 4. 11·106 0 0. 3 0. 34

**Period II.**

2. 98·108 1. 67·107 3. 02·108 8. 2·108 – 0. 03 – 0. 12 – 0. 03

**Period III**

2. 91·108 1. 3·106 5. 72·108 5. 81·105 0. 3 – 0. 005 0. 21

**Table 1.** Distribution of correlation coefficients during the 1 – III periods during the determined time

Days Period 129 – 182 Period 203 – 259 Period 307 – 385

Northern part R(Ea, Ed)=0. 02 R(Ea, Ed)=0. 52

**Table 2.** Distribution of correlation coefficients during the IV period during the determined time intervals.

Southern part R(Ea, Ed)=0. 68 R(Ea, Ed)=0. 24

(A, B)

Northern part(A) Southern part (B) Correlation coefficients

absorbed energy.

intervals.

**Figure 10.** Distribution of absorbed (1 – 2) and dissipated (3 – 4) energy of the whole mine field during the period III 24. 05. 2008 – 26. 07. 2009. The horizontal axes is time per days

**Figure 11.** Distribution of absorbed (1 – 2) and extracted (3 – 6) energy of the whole mine field during the period IV 28. 06. 2009 – 18. 07. 2010. The horizontal axes is time per days

The whole interval of research had been divided into four periods: from 03. 06. 2006 – 13. 01. 2007 (period I), from 14. 01. 2007 – 17. 05. 2008 (period II), from 24. 05. 2008 – 26. 07. 2009 (period III), from 28. 06. 2009 – 18. 07. 2010 (period IV). (Figures 8 – 11).

98 Fractal Analysis and Chaos in Geosciences

**Figure 10.** Distribution of absorbed (1 – 2) and dissipated (3 – 4) energy of the whole mine field during

**Figure 11.** Distribution of absorbed (1 – 2) and extracted (3 – 6) energy of the whole mine field during

the period IV 28. 06. 2009 – 18. 07. 2010. The horizontal axes is time per days

the period III 24. 05. 2008 – 26. 07. 2009. The horizontal axes is time per days

We would like to note such peculiarities for the fourth period (Figure 11): during the period 100 days, beginning from 100 days from the beginning of the analyzed period and finishing by 200 days, the explosions had been provided in the northern and in the southern part of the mine field approximately of equal intensity, however the energy of massif response in the southern part is significantly larger, than from the northern part. During next 50 days the explosions had been provided in the southern part, but the energy of the massif response in the northern part and in the southern part are approximately equal. During the period from 300 days to 400 days the explosions had been provided mainly in the northern part of the mine field, the distribution of the massif response energy in the northern part practically corresponds to the distribution of the absorbed energy.

**Table 1.** Distribution of correlation coefficients during the 1 – III periods during the determined time intervals.

**Table 2.** Distribution of correlation coefficients during the IV period during the determined time intervals.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 101

**Figure 13.** Geoelectrical sections along the profile horizon – 210, NWP, ort 3, frequency 10. 16 kHz. a) 9-th of

July 2010, b) 16-th of July 2010, c) 23-th of July 2010, d) 26-th of July 2010, *M*0

**Figure 12.** Phase diagrams: a) period 1, b) period 2, c) period 3, d) period 4., 1 – southern part, 2 – northern part.

2 – northern part.

**Figure 12.** Phase diagrams: a) period 1, b) period 2, c) period 3, d) period 4., 1 – southern part,

**Figure 13.** Geoelectrical sections along the profile horizon – 210, NWP, ort 3, frequency 10. 16 kHz. a) 9-th of July 2010, b) 16-th of July 2010, c) 23-th of July 2010, d) 26-th of July 2010, *M*0 -intensity of discrete zones .

Here, looking the table 1, we can see the changing of nonlinearity in time of the massif state by its active blast effects. The practical massif state situation cannot be described by a linear model. Phase diagrams reflect the dependence of massif energy dissipation velocity from the irregularity in time active man-caused influence. So during the first period 03. 06. 2006 – 13. 01. 2007 the influence of that irregularity leads to dynamical events for the southern part 27. 08 2006. и 17. 09. 2006., in the northern part 10. 09. 2006, however the intensity of the dynamical events during that period did not exceed, than 5 ·106 joules. During the second period 14. 01. 2007 – 17. 05. 2008 the influence of the irregularity significantly increases compared to the first period, especially in the northern part of the mine and there occurred some dynamical events on 13. 05, 12. 08 and on 21. 10. 2007 in the southern part the most powerful rock burst occurred with energy amplitude – 109 joules.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 103

**Figure 14.** The distribution of absorbed (Ep) and dissipated energy (Ev) by the first and second part of the southern part of the massif during the second period of observation time 14. 01. 2007 – 17. 05. 2008.

Intervals R(lgEp, lgEv1) R(lgEp, lgEv2) R(lgEv2, lgEv1)

207 – 336 – 0. 00881 0. 12652

105 – 207 0. 30794 0. 13077 0. 70973

0 – 200 0. 26314

**Table 3.** Coefficients of correlation R for processes between absorption energy (Ep) and dissipation

Here we divided the southern part of the mine's field on two parts: orts 13 – 16 – block (1) and 23 – 31 – block (2), horizons – 280, – 350 and used the data from the mines catalogue from 14-th January 2007 to 17-th of May 2008 (period 2).. From the results of the table 3 we see that the process of energy adsorption and dissipation in the researched blocks 1 and 2 as a rule is nonlinear, but the degree of nonlinearity changes in time. In the interval 23 – 42, which corresponds to the interval DT (200 – 342), the correlation coefficient between lg(Ep)(DT) and lg(Ev2)(DT) has the maximum value. That interval includes the process of

energy (Ev) in the massif for different time intervals (Ni),for the second period.

The horizontal ax is time per days.

1 – 21

21 – 60

23 – 42

24 – 41

12 – 24

1 – 23

Intervals R(lgEp, lgEv2)

0 – 179 – 0. 0062

179 – 450 0. 240683

200 – 342 0535269

The phase diagrams of the massif state during the two last periods show the identity of the massif state of the southern and northern parts during the last two years. The energy of responses and the velocity of its changing's has an identical character.

The successive cycle of induction active electromagnetic monitoring was provided in 2010 from 9-th of July to 26-th of July in the holes of the northern –west department and in some holes of the southern part of the mine's field. During the same time the man-caused works had been arranged in the block 4 – 5 of the northern part of the mine. The explosions had been achieved 04. 07 – Energy of the explosion 1. 7Е+ 06 joules, 11. 07 – Energy of the explosion 3. 50E+ 07 joules, 18. 07 – Energy of the explosion 1. 7Е+ 06 joules, 25. 07 – Energy of the explosion 1. 7Е+ 06 joules, 01. 08 – Energy of the explosion 1. 4Е+ 05 joules. The repeated electromagnetic observations had been achieved in the ort 3 NWD 9. 07, 16. 07, 23. 07 and 26. 07. (figure 10, (a-d)) The analyze of electromagnetic data during the 4 cycles of observation from 2007 to 2010 showed that the massif of the 3-d ort is more sensitive to the change of the stress-deformed state in the northern-west part of the mine field, which is caused by the influence outside it.

## **4. Geosynergetics approach for analyze of rock state – Theoretical and experimental results**

The research of rock burst hazard massif of the mine Tashtagol was arranged using the approaches of the theory of open dynamical systems [6, 8, 20]. We would like to search the criterions of changing the regimes of dissipation for the real rock massif, which are under heavy man-caused influence. We used the data from the seismic catalogue from June 2006 to June 2008, used the space-time coordinates of whole dynamical eventsresponses of the massif, which occurred during that period and also explosions, which are fixed by the seismic station using the energy parameter [20]. The phase portraits of the massif state of the northern and southern places had been developed in the coordinates Ev (t) and d(Ev(t))/dt, t-time in the parts of days, Ev-seismic energy dissipated by the massif in joules.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 103

102 Fractal Analysis and Chaos in Geosciences

caused by the influence outside it.

**experimental results**

in joules.

Here, looking the table 1, we can see the changing of nonlinearity in time of the massif state by its active blast effects. The practical massif state situation cannot be described by a linear model. Phase diagrams reflect the dependence of massif energy dissipation velocity from the irregularity in time active man-caused influence. So during the first period 03. 06. 2006 – 13. 01. 2007 the influence of that irregularity leads to dynamical events for the southern part 27. 08 2006. и 17. 09. 2006., in the northern part 10. 09. 2006, however the intensity of the dynamical events during that period did not exceed, than 5 ·106 joules. During the second period 14. 01. 2007 – 17. 05. 2008 the influence of the irregularity significantly increases compared to the first period, especially in the northern part of the mine and there occurred some dynamical events on 13. 05, 12. 08 and on 21. 10. 2007 in the southern part the most

The phase diagrams of the massif state during the two last periods show the identity of the massif state of the southern and northern parts during the last two years. The energy of

The successive cycle of induction active electromagnetic monitoring was provided in 2010 from 9-th of July to 26-th of July in the holes of the northern –west department and in some holes of the southern part of the mine's field. During the same time the man-caused works had been arranged in the block 4 – 5 of the northern part of the mine. The explosions had been achieved 04. 07 – Energy of the explosion 1. 7Е+ 06 joules, 11. 07 – Energy of the explosion 3. 50E+ 07 joules, 18. 07 – Energy of the explosion 1. 7Е+ 06 joules, 25. 07 – Energy of the explosion 1. 7Е+ 06 joules, 01. 08 – Energy of the explosion 1. 4Е+ 05 joules. The repeated electromagnetic observations had been achieved in the ort 3 NWD 9. 07, 16. 07, 23. 07 and 26. 07. (figure 10, (a-d)) The analyze of electromagnetic data during the 4 cycles of observation from 2007 to 2010 showed that the massif of the 3-d ort is more sensitive to the change of the stress-deformed state in the northern-west part of the mine field, which is

**4. Geosynergetics approach for analyze of rock state – Theoretical and**

The research of rock burst hazard massif of the mine Tashtagol was arranged using the approaches of the theory of open dynamical systems [6, 8, 20]. We would like to search the criterions of changing the regimes of dissipation for the real rock massif, which are under heavy man-caused influence. We used the data from the seismic catalogue from June 2006 to June 2008, used the space-time coordinates of whole dynamical eventsresponses of the massif, which occurred during that period and also explosions, which are fixed by the seismic station using the energy parameter [20]. The phase portraits of the massif state of the northern and southern places had been developed in the coordinates Ev (t) and d(Ev(t))/dt, t-time in the parts of days, Ev-seismic energy dissipated by the massif

powerful rock burst occurred with energy amplitude – 109 joules.

responses and the velocity of its changing's has an identical character.

**Figure 14.** The distribution of absorbed (Ep) and dissipated energy (Ev) by the first and second part of the southern part of the massif during the second period of observation time 14. 01. 2007 – 17. 05. 2008. The horizontal ax is time per days.

**Table 3.** Coefficients of correlation R for processes between absorption energy (Ep) and dissipation energy (Ev) in the massif for different time intervals (Ni),for the second period.

Here we divided the southern part of the mine's field on two parts: orts 13 – 16 – block (1) and 23 – 31 – block (2), horizons – 280, – 350 and used the data from the mines catalogue from 14-th January 2007 to 17-th of May 2008 (period 2).. From the results of the table 3 we see that the process of energy adsorption and dissipation in the researched blocks 1 and 2 as a rule is nonlinear, but the degree of nonlinearity changes in time. In the interval 23 – 42, which corresponds to the interval DT (200 – 342), the correlation coefficient between lg(Ep)(DT) and lg(Ev2)(DT) has the maximum value. That interval includes the process of preparation of the resonance release energy by the massif as a rock burst of the 9-th class. From the other side for the interval 12 – 42 we see the changing of the type of the correlation dependence between the functions R(lg(Ev2), lg(Ev1)): for the interval 12 – 24 the correlation coefficient corresponds to a linear correlation function, which reflects the elastic interaction between two blocks (orts 15 – 16) and (orts 27 – 29). By the way the relation between lg (Ep) and lg (Ev2) and lg(Ev1) is practically nonlinear. For the interval 24 – 41 the type of correlation function R(lg(Ep2), lg(Ep1)) changes nonlinear and the dependence between lg(Ep) and lg(Ev1) is practically absent (table 3). However before 48 days a rock burst of 6.4 th class occurs in the block 1 after the explosion of 5.2 classes in the block 2. Can we think that this rock burst was a foreshock for the rock shock of the 9-th class in the block 2? On that question we can answer only after providing detailed every day observations of electromagnetic active induction and deformation monitoring in the whole space of the 1-st and second blocks.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 105

further movement of phase points lead to return to the same attracting area. These movements may take place in resonance or spontaneous cases. The received practical results are significant because firstly we could find the coincidences with the mathematical theory between open systems and results of experiments in natural medium with very complicated structure. On this base we developed a new processing method for the seismological information which can be used real-time for estimation of the disaster

On the base of the constructed algorithm for calculation of the distribution of seismic waves into a medium with hierarchic inclusions located in an arbitrary layer of a horizontal layered elastic medium we can compute the stress components on each hierarchic level. This information we use to estimate the medium state analyzing the hierarchic structure and its changing. From the other side as higher the degree of hierarchic structure as larger becomes the degree of space nonlinearity of seismic field components distribution. This feature should be taken into account by interpretation to minimize linearization negative effects. From the received theoretical results we did a conclusion, that as higher is rank of hierarchy of the medium as less the similarity between seismic and electromagnetic results, and the obtained information has a independent sense, which underlines the complexity of the

The analysis of experimental seismological and electromagnetic information showed the common additional information on different space-time scale levels of the state of rock massif, which is under energetic influence of mining explosions. It has revealed the change of nonlinearity in time of the massif state. The description of the massif behavior in the frame of the linear dynamical model does not correspond to the real practical situation. As it follows from the received results, the changing of the massif state: decreasing and increasing of its activization, does not depend on the space location of the explosion and on delay in time. It is needed to determine that delay function of its activization to be able to forecast the massif behavior. For that it is needed to continue obtaining and analyzing the complex information from passive and active seismic and electromagnetic detailed monitoring. For quantitative research of the behavior of different types of nonlinear dynamical systems we shall use the asymptotical methods for nonlinear mechanics, developed by N. M. Krilov, N. N. Bogolubov and Yu. A. Mitropolsky [bibliography in [23]]. However if we want to use the data from the seismic catalogue, we must convert energetic characteristics to forces and displacements. For that we need additional information about deformations, which occur in the massif by the explosions influence. The mathematical method, which had been developed by academician N. N. Bogolubov, allows us to get on with quantitative description of the causes of self excitation of the nonlinear mechanical system and the occurrence of the space-time local resonance in the system as a response to the outer influence. Using of common theoretical approaches [8, 23] and complex data: seismic catalogue data, deformation and induction electromagnetic monitoring data [20, 22] will allow us to formulate and solve the problem of forecasting the critical state of activated local

danger degree changing in mine massif.

researched medium.

place in rock massif.

There is a deep principal difference between the mechanics of linear and nonlinear oscillations, which still persists by research weak nonlinear oscillations, which can be described by differential equations, which differ from linear equations only by presence of small terms, which begin to influence especially on the time intervals larger, than the period of oscillations. In the system there can be energy sources and absorbers, which produce and absorb very small work for one period of oscillations, but by prolonged actions their produced effect can be summarized and provide a sufficient influence on the oscillatory processing: decay, increase, stability. The small nonlinear terms can provide cumulative influence and damage the superposition principle, apart harmonics begin to influence on each other, and therefore the individual research of the behavior each harmonic oscillation apart cannot be done. The continuous waves can practically exist only in the case, if the system contains an energy source, which can compensate the energy decay, which appears as a result of existence of dissipative forces. Such source plays a role of negative friction. The oscillations which occur owing to a source, which influence has not constant period, are known as auto oscillations and any auto oscillatory system is described by a non linear differential equation.

Relaxation oscillations are widely distributed in the nature, for them oscillatory process has two stages: slow energy accumulation by the system and then energy relaxation, which take place almost immediately after the moment when potential threshold is over reached for that system.

## **5. Conclusions**

One of the mathematical ideas about the common causes of chaotization and stochastization of dynamical system movements are its losses of stability and exponential recession of near located phase trajectories together with its common boundedness and its common compression [8]. The mathematical result coincides as a whole with the practical result: in the phase area the smaller attracting phase trajectories area exists where may occur an exponential recession of them, then the movement character changes and the further movement of phase points lead to return to the same attracting area. These movements may take place in resonance or spontaneous cases. The received practical results are significant because firstly we could find the coincidences with the mathematical theory between open systems and results of experiments in natural medium with very complicated structure. On this base we developed a new processing method for the seismological information which can be used real-time for estimation of the disaster danger degree changing in mine massif.

104 Fractal Analysis and Chaos in Geosciences

and second blocks.

differential equation.

that system.

**5. Conclusions**

preparation of the resonance release energy by the massif as a rock burst of the 9-th class. From the other side for the interval 12 – 42 we see the changing of the type of the correlation dependence between the functions R(lg(Ev2), lg(Ev1)): for the interval 12 – 24 the correlation coefficient corresponds to a linear correlation function, which reflects the elastic interaction between two blocks (orts 15 – 16) and (orts 27 – 29). By the way the relation between lg (Ep) and lg (Ev2) and lg(Ev1) is practically nonlinear. For the interval 24 – 41 the type of correlation function R(lg(Ep2), lg(Ep1)) changes nonlinear and the dependence between lg(Ep) and lg(Ev1) is practically absent (table 3). However before 48 days a rock burst of 6.4 th class occurs in the block 1 after the explosion of 5.2 classes in the block 2. Can we think that this rock burst was a foreshock for the rock shock of the 9-th class in the block 2? On that question we can answer only after providing detailed every day observations of electromagnetic active induction and deformation monitoring in the whole space of the 1-st

There is a deep principal difference between the mechanics of linear and nonlinear oscillations, which still persists by research weak nonlinear oscillations, which can be described by differential equations, which differ from linear equations only by presence of small terms, which begin to influence especially on the time intervals larger, than the period of oscillations. In the system there can be energy sources and absorbers, which produce and absorb very small work for one period of oscillations, but by prolonged actions their produced effect can be summarized and provide a sufficient influence on the oscillatory processing: decay, increase, stability. The small nonlinear terms can provide cumulative influence and damage the superposition principle, apart harmonics begin to influence on each other, and therefore the individual research of the behavior each harmonic oscillation apart cannot be done. The continuous waves can practically exist only in the case, if the system contains an energy source, which can compensate the energy decay, which appears as a result of existence of dissipative forces. Such source plays a role of negative friction. The oscillations which occur owing to a source, which influence has not constant period, are known as auto oscillations and any auto oscillatory system is described by a non linear

Relaxation oscillations are widely distributed in the nature, for them oscillatory process has two stages: slow energy accumulation by the system and then energy relaxation, which take place almost immediately after the moment when potential threshold is over reached for

One of the mathematical ideas about the common causes of chaotization and stochastization of dynamical system movements are its losses of stability and exponential recession of near located phase trajectories together with its common boundedness and its common compression [8]. The mathematical result coincides as a whole with the practical result: in the phase area the smaller attracting phase trajectories area exists where may occur an exponential recession of them, then the movement character changes and the On the base of the constructed algorithm for calculation of the distribution of seismic waves into a medium with hierarchic inclusions located in an arbitrary layer of a horizontal layered elastic medium we can compute the stress components on each hierarchic level. This information we use to estimate the medium state analyzing the hierarchic structure and its changing. From the other side as higher the degree of hierarchic structure as larger becomes the degree of space nonlinearity of seismic field components distribution. This feature should be taken into account by interpretation to minimize linearization negative effects. From the received theoretical results we did a conclusion, that as higher is rank of hierarchy of the medium as less the similarity between seismic and electromagnetic results, and the obtained information has a independent sense, which underlines the complexity of the researched medium.

The analysis of experimental seismological and electromagnetic information showed the common additional information on different space-time scale levels of the state of rock massif, which is under energetic influence of mining explosions. It has revealed the change of nonlinearity in time of the massif state. The description of the massif behavior in the frame of the linear dynamical model does not correspond to the real practical situation. As it follows from the received results, the changing of the massif state: decreasing and increasing of its activization, does not depend on the space location of the explosion and on delay in time. It is needed to determine that delay function of its activization to be able to forecast the massif behavior. For that it is needed to continue obtaining and analyzing the complex information from passive and active seismic and electromagnetic detailed monitoring. For quantitative research of the behavior of different types of nonlinear dynamical systems we shall use the asymptotical methods for nonlinear mechanics, developed by N. M. Krilov, N. N. Bogolubov and Yu. A. Mitropolsky [bibliography in [23]]. However if we want to use the data from the seismic catalogue, we must convert energetic characteristics to forces and displacements. For that we need additional information about deformations, which occur in the massif by the explosions influence. The mathematical method, which had been developed by academician N. N. Bogolubov, allows us to get on with quantitative description of the causes of self excitation of the nonlinear mechanical system and the occurrence of the space-time local resonance in the system as a response to the outer influence. Using of common theoretical approaches [8, 23] and complex data: seismic catalogue data, deformation and induction electromagnetic monitoring data [20, 22] will allow us to formulate and solve the problem of forecasting the critical state of activated local place in rock massif.

That work is fulfilled by support of RFBR grant 10 – 05 – 00013a

## **Author details**

O.A. Khachay\* *Institute of Geophysics UB RAS, Yekaterinburg, Russia*

A.Yu. Khachay and O.Yu. Khachay *Ural Federal University, Institute of Mathematics and Computer Sciences, Yekaterinburg, Russia*

## **6. References**

[1] Physical mesomechanics and computer construction of materials. Novosibirsk: Nauka, SIFR; 1995. V. 1.

Dynamical Model for Evolution of Rock Massive State as a Response on a Changing of Stress-Deformed State 107

[11] Dmitriev V. I. The diffraction of the plane electromagnetic field on cylindrical bodies, located in layered medium. In: Voevodin V. V. Methods of calculation and

[12] Kupradze V. D. The boundary problems of the theory of waves and integral equations.

[13] Khachay A. Yu. Algorithm of direct problem solution of dynamical seismic research by excitation with a vertical point force source, located in an arbitrary layer of n-layered isotropic elastic medium. In: Mazurov Vl. D., A. I. Smirnov (eds) Informatics and

[14] Hachay O. A., Novgorodova E. N., Khachay A. Yu. Research of resolution of planshet electromagnetic method for active mapping and monitoring of heterogeneous

[15] Hachay O. A., Hinkina T. A., Bodin V. V. Preconditions of seismic and electromagnetic monitoring of nonstationary medium. Russian geophysical journal 2000; 17 – 18 83 – 89. [16] Hachay O. A., Hinkina T. A., Bodin V. V. Research of the criterion of similarity for seismic and electromagnetic research in frequency-geometrical variant. In: Holshevnikov K. V., Shustov B. M. Zacharova P. E. (eds) Astronomical and geodetic

[17] Hachay O. A., Druginin V. S., Karetin Yu. S., Bodin V. V., Novgorodova E. N., Khachay A. Yu., Zacharov I. B., Groznih M. V. The use of complex planshet seismic and electromagnetic method for solution of mapping subsurface heterogeneities problems. In: Solodilov L. N., Kostuchenko S. L., Yasuleyitch N. N. (eds.): Geophysics of the 21 – st century: proceedings of the third geophysical symposium, 22 – 24 February 2001,.

[18] Hachay. O. A. Mathematical modeling and interpretation of alternative electromagnetic field for heterogeneous crust and upper mantle of the Earth. Prof. thesis. Yekaterinburg:

[19] Hachay O. A., Khachay O. Yu. Simulation of seismic and electromagnetic field in hierarchic heterogeneous structures. In: Martyshko P. S. (Ed): Geophysical research of the Ural and adjusted regions: proceedings of the international conference, 4 – 8

[21] Hachay O. A., Khachay A. Yu. Construction of dynamical Model for Evolution of Rock Massive State as a Response on a Changing of stress-deformed State. Geophysical

[22] Hachay O. A., Khachay O. Yu, Klimko V. K., Shipeev O. Yu. The reflection of synergetic features in the response of geological medium on outer force actions. In: Jinghong F., Junqian Z., Haibo C., Zhaohui J.: Advances in heterogeneous Material Mechanics:

February 2008, Yekaterinburg, Russia. Yekaterinburg: OOO"IRA UTK", 2008. [20] Hachay O. A., Khachay O. Yu., Klimko V. K., Shipeev O. V. Reflection of synergetic features of rock massif state under man-made influence in the data of seismological

catalogue. Mining information and analytical bulletin 2010; 6 259 – 271.

programming in layered medium. Moscow: MSU, 1965. V. 3. p. 307 – 316.

mathematical modeling. Yekaterinburg: USU; 2006. p. 279 – 310.

geoelektrics medium. Fizika Zemli 2003; 1 30 – 41.

research. Yekaterinburg: USU; 2001. p. 30 – 35.

Moscow, Russia. Moscow: Nauchniy Mir; 2001.

Research abstracts 2010; 12. EGU2010 – 2662.

IGF UB RAS; 1994.

Moscow-Leningrad: TTL. 1950,

\* Corresponding Author

[11] Dmitriev V. I. The diffraction of the plane electromagnetic field on cylindrical bodies, located in layered medium. In: Voevodin V. V. Methods of calculation and programming in layered medium. Moscow: MSU, 1965. V. 3. p. 307 – 316.

106 Fractal Analysis and Chaos in Geosciences

A.Yu. Khachay and O.Yu. Khachay

**Author details**

*Ural Federal University,*

SIFR; 1995. V. 1.

2006; 5 109 – 115.

**6. References**

2003.

 \*

Corresponding Author

"LIBROKOM"; 2009.

O.A. Khachay\*

That work is fulfilled by support of RFBR grant 10 – 05 – 00013a

*Institute of Mathematics and Computer Sciences, Yekaterinburg, Russia*

[1] Physical mesomechanics and computer construction of materials. Novosibirsk: Nauka,

[2] Hachay O. A. The problem of the research of redistribution of stress and phase states of massif between high man-made influences. Mining information and analytic bulletin

[3] Hachay O. A., Khachay O. Yu. Theoretical approaches for system of geophysical state control validation of geological medium by man-made influence. Mining information

[4] Hachay O. A., Khachay O. Yu. Results of electromagnetic and seismic monitoring of the state of rock massif by use the approach of the open dynamical systems. Geophysical

[5] Hachay O. A. Synergetic events in geological medium and nonlinear features of wave

[6] Chulichkov A. I. Mathematical models of nonlinear dynamics. Moscow: Phismatlit;

[8] Naimark Yu. I., Landa P. S. Stochastic and chaotic oscillations. Moscow; Knigniy dom

[9] Hachay O. A. Geosynergetics: theory, methods, experiments. In: Spichak V. V. (ed.)

[10] Hachay O. A., Khachay A. Yu. Theory of integrating seismic and electromagnetic active methods for mapping and monitoring the state of 2D heterogeneities with hierarchic structure in a N-layered medium. In: Shkuratnik V. L., Fokin A. V., Mironov A. V., Didenkulov I. N. (eds) Physical acoustics, nonlinear acoustics, distribution and wave diffraction, geoacoustics : proceedings of the XX-th session of the Russian acoustic

Electromagnetic research of the Earth. Moscow: Krasand; 2009. p. 138 – 153.

society 27 – 31 October 2008, Moscow. Russia. VI. Moscow: GEOS; 2008.

propagation. Geophysical research abstracts 2009; 11 EGU2009 – 3684.

[7] Malineckiy G. G. Mathematical base of synergetics. Moscow; LKI; 2007.

*Institute of Geophysics UB RAS, Yekaterinburg, Russia*

and analytic bulletin 2008; 1 161 – 169.

research abstracts 2009; 11 EGU2009 – 137.

proceedings of the third international conference on heterogeneous material mechanics 22 – 26 May 2011, Shanghai, China. Lancaster, USA: DEStech Publications, Inc. 2011. [23] Bogolubov N. N. Mathematics and nonlinear mechanics. Moscow: Nauka; 2005.

**Chapter 6**

© 2012 Biswas et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Biswas et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Application of Multifractal and Joint Multifractal**

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

 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

 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

**Analysis in Examining Soil Spatial Variation:**

Asim Biswas, Hamish P. Cresswell and Bing C. Si

Additional information is available at the end of the chapter

**A Review**

http://dx.doi.org/10.5772/51437

**1. Introduction**

soil spatial variation.

spatial autocorrelation or similarity.

variability is also scale dependent.

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

Asim Biswas, Hamish P. Cresswell and Bing C. Si

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/51437

## **1. Introduction**

108 Fractal Analysis and Chaos in Geosciences

proceedings of the third international conference on heterogeneous material mechanics 22 – 26 May 2011, Shanghai, China. Lancaster, USA: DEStech Publications, Inc. 2011.

[23] Bogolubov N. N. Mathematics and nonlinear mechanics. Moscow: Nauka; 2005.

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.

© 2012 Biswas et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Biswas et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

 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.

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

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

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 *Z* as *q z* , 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

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

stationarity assumption of geostatistics.

and simulate spatially-variable soil properties.

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.

110 Fractal Analysis and Chaos in Geosciences

another.

locations [2, 22].

behaviour is known as the periodicity.

adding the individual observational results together.

 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

 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

 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

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

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 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 *Z* as *q z* , 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.

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

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

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

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),

( ) ( ) *q q Zx x*

where the symbol " " indicates proportionality, *Z* is the spatial series, and *x* is the lag distance. If the plot of ߬ሺݍሻ vs. *q* [or ߬ሺݍሻ curve; Fig. 1a] has a single slope (i.e., a straight line), then the series is a simple scaling (monofractal) type. If the ߬ሺݍሻ 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 ߬ሺݍሻ 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

the data are scaling in that range—i.e.,

(1)

( ) *q* , is

examining the soil spatial variation.

defined as [46]:

follows a power law equation of the form *h h*

there is certain degree of statistical scale-invariance [52].

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 science to 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 therefore be 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 science within 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 as a transect) as the mass (or value) distribution of a statistical

brief Conclusions Section.

**2.1. Multifractal analysis**

**2. Methods**

over a range of spatial and temporal scales.

quantitative framework for integrating soil biological, chemical and physical phenomena

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 science to 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 therefore be 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 science within a

Multifractal analysis can be used to characterize the scaling property of a soil variable measured in a direction (such as a 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 *h h* 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), ( ) *q* , is defined as [46]:

$$
\langle \left[ \Delta Z(\mathbf{x}) \right]^q \rangle \propto \mathbf{x}^{\mathbf{r}(q)} \tag{1}
$$

where the symbol " " indicates proportionality, *Z* is the spatial series, and *x* is the lag distance. If the plot of ߬ሺݍሻ vs. *q* [or ߬ሺݍሻ curve; Fig. 1a] has a single slope (i.e., a straight line), then the series is a simple scaling (monofractal) type. If the ߬ሺݍሻ 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 ߬ሺݍሻ 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 ( ) *q* curve and the UM model using the chi-square test. Statistically significant difference between the curves indicates multifractal scaling and nonsignificant difference indicates monofractal scaling. The degree of monofractal/multifractal can also be examined by calculating the deviation of the ( ) *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 ( ) *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 ( ) *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( ) *q* curve and thus multifractal behavior [47, 53].

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

( ) *q*

organic carbon, etc.) in a given segment relative to the whole support length. The

used widely in multifractal analysis.

and the local scaling indices, ߙ, are given by

calculated to the limit as

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

The multifractal spectrum, *f q* (Fig. 1b), which is the fractal dimension of the subsets of segments of size ߝ units with a coarse Hölder exponent (local scaling indices) of ߙ if

1

( ) log ( , )log ( , ) lim *i i i f q q q*

1

( ) log ( , )log ( ) lim *i i i q q p*

*f q qq q* ( ( )) ( ) ( )

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

For many practical applications indicator parameters are selected and used, in addition to the multifractal spectrum (*f*(*q*) vs. ߙ)*q*)), to describe the scaling property and variability of a process. Generalized dimensions, *Dq* , (Fig. 1c) are often used to provide indicator

0

The *Dq* value at *q* = 0, *D*0, is called the capacity dimension or the box counting dimension of the geometric support of the measure. The value at *q* = 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*1 close to unity indicates the evenness of measures over the sets of a given cell size,

log ( ) 1 lim 1 log( )

*i i*

*p*

(3)

(5)

(6)

(4)

0 , can be expressed as [36, 50],

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

0

0

*L*

largest value is the capacity dimension that is obtained at *q* = 0.

parameters. The *Dq* of a multifractal measure is calculated as

*q*

*q*

*D*

*L*

**Figure 1.** A typical of a) mass exponent graph (or ( ) *q* curve), b) multifractal spectrum (or *f q* 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 *qth* order normalized probability measures of a variable (also known as the partition function), (,) *q* , vary with the scale size, in a manner similar to Eq. (1) [36, 51], i.e.,

$$\mu\_i(q,\varepsilon) = \frac{\left[\,\_p\mu\_i(\varepsilon)\right]^q}{\sum\_i \left[\,\_p\mu\_i(\varepsilon)\right]^q} \propto \left(\varepsilon / L\right)^{\varepsilon(q)}\tag{2}$$

where ( ) *i p* is the probability of a measure in the *th i* segment of size units, calculated by dividing the value of the variable in the segment by the whole support length of *L* units. In simpler terms, ( ) *i p* 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 ( ) *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 q* (Fig. 1b), which is the fractal dimension of the subsets of segments of size ߝ units with a coarse Hölder exponent (local scaling indices) of ߙ if calculated to the limit as 0 , can be expressed as [36, 50],

$$f(q) = \lim\_{\varepsilon \to 0} \left( \log \left( \frac{\varepsilon}{L} \right) \right)^{-1} \sum\_{i} \mu\_i(q, \varepsilon) \log \mu\_i(q, \varepsilon) \tag{3}$$

and the local scaling indices, ߙ, are given by

114 Fractal Analysis and Chaos in Geosciences

goodness-of-fit between the

the slopes indicates nonlinearity in the

**Figure 1.** A typical of a) mass exponent graph (or

behavior in any property.

Eq. (1) [36, 51], i.e.,

where ( ) *i p*

simpler terms, ( ) *i p*

known as the partition function),

line fitted to the

can also be examined by calculating the deviation of the

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

Statistically significant difference between the curves indicates multifractal scaling and nonsignificant difference indicates monofractal scaling. The degree of monofractal/multifractal

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

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

Research has indicated that the *qth* order normalized probability measures of a variable (also

( ) ( ) (,) (/) ( )

dividing the value of the variable in the segment by the whole support length of *L* units. In

*p q L p*

*i q i i*

*q i q*

, vary with the scale size,

measures the concentration of a variable of interest (e.g., clay content,

*i* segment of size

and c) generalized dimension graph (or *Dq* curve) used to identify the monofractal and multifractal

(,) *q*

is the probability of a measure in the *th*

UM model. Significant difference in the slope indicates multifractal behavior. The

( ) *q* curve and the UM model using the chi-square test.

( ) *q* curve and thus multifractal behavior [47, 53].

( ) *q* curve), b) multifractal spectrum (or *f q* curve)

in a manner similar to

units, calculated by

(2)

( ) *q* curve from the UM model [47,

( ) *q* can

( ) *q* curve (also referred as single fit) can be compared to the slopes of the

$$\alpha(q) = \lim\_{\varepsilon \to 0} \left( \log \left( \frac{\varepsilon}{L} \right) \right)^{-1} \sum\_{i} \mu\_i(q, \varepsilon) \log p\_i(\varepsilon) \tag{4}$$

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

$$f(\alpha(q)) = q\alpha(q) - \tau(q) \tag{5}$$

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. ߙ)*q*)), to describe the scaling property and variability of a process. Generalized dimensions, *Dq* , (Fig. 1c) are often used to provide indicator parameters. The *Dq* of a multifractal measure is calculated as

$$D\_q = \frac{1}{q - 1} \lim\_{\varepsilon \to 0} \frac{\log \sum\_i p\_i(\varepsilon)}{\log(\varepsilon)}\tag{6}$$

The *Dq* value at *q* = 0, *D*0, is called the capacity dimension or the box counting dimension of the geometric support of the measure. The value at *q* = 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*1 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*1 and *D*2 become similar to the *D*0, the capacity dimension. The value of *D*0 = *D*1 = *D*2 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*0>*D*1>*D*2. The *D*1/*D*0 value can also be used to describe the heterogeneity in the distribution [56]. A *D*1/*D*0 value equal to 1 indicates exact monoscaling of the distribution, which means that all fractions take equal values at different scales.

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

() ()

*p r*

(,) *q t* and

*q t j j*

*i i*

*i i*

*i i*

(8)

(9)

(,) *q t* and

(10)

*q t i i*

  , where *α* and *β* are

, respectively. The partition

(7)

(,) *q t* respectively, are calculated as

(,) *q t* are the mean local

 and *Ri* ,

 and *Ri*

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

and for the second variable is *R L i*

( )

() () (,, )

*p r q t*

1

The local scaling indices (coarse Hölder exponents) with respect to two probability measures

( ) 1 1 ( , ) log( ( )) ( , , ) log( ( )) *N*

*i qt N qt p*

( ) 1 1 ( , ) log( ( )) ( , , ) log( ( )) *N*

*i qt N qt r*

( ) 1 1 ( , ) log( ( )) ( , , ) log( ( , , )) *N*

*i f N qt qt*

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 ݂ሺߙǡ ߚሻ is attained, which is the dimension if all the segments contain the same concentration of mass. Different pairs of ߙ and ߚ 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

*j*

function (the normalized *μ* measures) for the joint distribution of *i P*

*i N*

) of the set on which

weighted by the real numbers *q* and *t* can be calculated as [47, 50, 51];

, which are represented by

> >

exponents of both measures is given by

**2.2. Joint multifractal analysis**

the local singularity strength corresponding to *i P*

variable is *i P L*

( ) *i P*

follows:

 and ( ) *Ri*

The dimension (i.e., *f*(,)

#### *2.1.1. Multifractal analysis steps*

	- a. 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 ( ,) ( ) *i m i p M* where *i* is the index for *i* th segment when the whole transect is divided into *n* segments with a unit length = *i n* 1, , ; *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.

b. 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).

$$\text{a. } \quad \text{Calculate } \, F f(q\_{\prime} \varepsilon) = \sum\_{i} \mu\_{i}(q\_{\prime} \varepsilon) \log \mu\_{i}(q\_{\prime} \varepsilon).$$

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

#### **2.2. Joint multifractal analysis**

116 Fractal Analysis and Chaos in Geosciences

different scales.

*2.1.1. Multifractal analysis steps*

dataset with regular interval between samples.

divide by the total of all segments ( ,) ( )

length *ε*; *M* is the total of all points along the transect.

*i Fq q p*

3. Calculate *τ*(*q*) as the intercept of linear regression of log ( , ) *i*

6. Calculate *Dq* as the intercept of the linear regression of *FD q*(,)

7. Plot *τ*(*q*) as a function *q*, *f*(*q*) as a function of *α*(*q*), and *Dq* as a function of *q*.

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

4. Calculate *f*(*q*) as the intercept of linear regression of *Ff q*(,)

5. Calculate *α*(*q*) as the intercept of linear regression of *F q*

a. Calculate ( , ) ( , )log ( , ) *i i i Ff q q q*

b. Calculate ( , ) ( , )log ( ) *i i*

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*1 and *D*2 become similar to the *D*0, the capacity dimension. The value of *D*0 = *D*1 = *D*2 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*0>*D*1>*D*2. The *D*1/*D*0 value can also be used to describe the heterogeneity in the distribution [56]. A *D*1/*D*0 value equal to 1 indicates exact monoscaling of the distribution, which means that all fractions take equal values at

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

a. 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

*i*

b. For certain *q* values or statistical moments (e.g., -20 to 20, which are selected based

*i n* 1, , ; *m*(*ε*, *i*) is the sum of measurements at all points in the *i*

on the nature of the data to be analyzed), calculate *μi*(*q, ε*) using Eq. (2).

2. 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.

segment when the whole transect is divided into *n* segments with a unit length =

*m i p M*

 *q* vs. log( ).

(,)

 vs. log( ).

> vs. log( ).

> > vs. log( ).

where *i* is the index for *i*

th

th segment of

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 *i P L* and for the second variable is *R L i* , where *α* and *β* are the local singularity strength corresponding to *i P* and *Ri* , respectively. The partition function (the normalized *μ* measures) for the joint distribution of *i P* and *Ri* , weighted by the real numbers *q* and *t* can be calculated as [47, 50, 51];

$$\mu\_i(q, t, \varepsilon) = \frac{p\_i(\varepsilon)^q \cdot r\_i(\varepsilon)^t}{\sum\_{j=1}^{N(\varepsilon)} \left[ p\_j(\varepsilon)^q \cdot r\_j(\varepsilon)^t \right]} \tag{7}$$

The local scaling indices (coarse Hölder exponents) with respect to two probability measures ( ) *i P* and ( ) *Ri* , which are represented by (,) *q t* and (,) *q t* respectively, are calculated as follows:

$$\alpha(q, t) = -\left[\log(N(\varepsilon))\right]^{-1} \sum\_{i=1}^{N(\varepsilon)} \left[\mu\_i(q, t, \varepsilon) \cdot \log(p\_i(\varepsilon))\right] \tag{8}$$

$$\mathcal{J}(q,t) = -\left[\log(\mathcal{N}(\varepsilon))\right]^{-1} \sum\_{i=1}^{N(\varepsilon)} \left[\mu\_i(q,t,\varepsilon) \cdot \log(r\_i(\varepsilon))\right] \tag{9}$$

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

$$f(\boldsymbol{\alpha}, \boldsymbol{\beta}) = -\left[\log(\mathcal{N}(\boldsymbol{\varepsilon}))\right]^{-1} \sum\_{i=1}^{N(\boldsymbol{\varepsilon})} \left[\mu\_i(\boldsymbol{q}, t, \boldsymbol{\varepsilon}) \cdot \log(\mu\_i(\boldsymbol{q}, t, \boldsymbol{\varepsilon})) \right] \tag{10}$$

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 ݂ሺߙǡ ߚሻ is attained, which is the dimension if all the segments contain the same concentration of mass. Different pairs of ߙ and ߚ 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 ݂ሺߙǡ ߚሻ represents the frequency of the occurrence of a certain value of ߙ and a certain value of ߚ, high values of ݂ሺߙǡ ߚሻ signifies a strong association between the values of ߙ and the values of ߚ. By perturbing *q* and *t*, we can examine the association of similar values (highs vs. highs or lows vs. lows) of ߙ and ߚ as well as dissimilar values (highs vs. lows) of ߙ and ߚ. Generally, a contour graph is used to represent the joint dimensions, ݂ሺߙǡ ߚሻ 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 analysis steps*

	- a. 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 ( ,) ( ) *Y i Y m i p M* where *i* is the index for *i* th segment when the

whole transect is divided into *n* segments with a unit length = *i n* 1, , ; Similarly, we can calculate the probability measure for the variable *Z*: ( ) ( ) *Z m i r* .

*i*

*Z*

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

 *q t* vs. log( )

(,, )

(,, )

 vs. log( )

 vs. log( )

 vs. log( )

(Eq. (8)).

(Eq. (9)).

(Eq. (10)).

(Eq. (7)).

d. Calculate ( , , ) ( , , )log ( ) *i i*

e. Calculate ( , , ) ( , , )log ( ) *i i*

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

*i F qt qt p*

*i F qt qt r*

3. Calculate *f*(*q*, *t*) as the intercept of linear regression of *Ff q t* (,, )

4. Calculate *α*(*q*, *t*) as the intercept of linear regression of *F qt*

5. Calculate *β*(*q*, *t*) as the intercept of linear regression of *F qt*

**2.3. Comments on multifractal and joint multifractal analysis**

a transect. However, the support *ε* becomes an area or volume.

moments and how the higher moments change with scale.

may diverge and the process of division needs to be stopped.

can be analyzed by multifractal analysis.

properties on various soil processes.

**science**

2. Calculate *τ*(*q, t*) as the intercept of linear regression of log ( , , )

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

Multifractal analysis for two-dimensional or three dimensional fields is the same as for

 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

Multifractal analysis does not require regularly spaced samples. Any sampling scheme

 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

 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

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

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

*M*

d. Calculate ( , , ) ( , , )log ( ) *i i i F qt qt p*

118 Fractal Analysis and Chaos in Geosciences

vertical and horizontal axis (Fig. 2).

correlation and b) weak correlation.

*2.2.1. Joint multifractal analysis steps*

or a rectangle from an area.

using Eq. (7).

of all segments ( ,) ( ) *Y i*

similar moment orders. Because ݂ሺߙǡ ߚሻ represents the frequency of the occurrence of a certain value of ߙ and a certain value of ߚ, high values of ݂ሺߙǡ ߚሻ signifies a strong association between the values of ߙ and the values of ߚ. By perturbing *q* and *t*, we can examine the association of similar values (highs vs. highs or lows vs. lows) of ߙ and ߚ as well as dissimilar values (highs vs. lows) of ߙ and ߚ. Generally, a contour graph is used to represent the joint dimensions, ݂ሺߙǡ ߚሻ 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

**Figure 2.** A typical of multifractal spectra for joint distribution of two variables with a) strong

*Y m i*

whole transect is divided into *n* segments with a unit length =

*p M*

c. Calculate ( , , ) ( , , )log ( , , ) *i i i Ff q t q t q t*

1. Calculate the probability measure (*p*) for variable *Y* from a linear distance for a transect

a. 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

Similarly, we can calculate the probability measure for the variable *Z*: ( ) ( ) *Z*

b. For certain *q* values (e.g., -20 to 20) and t values (e.g., -20 to 20), calculate *μi*(*q, t, ε*)

where *i* is the index for *i*

th segment when the

*i n* 1, , ;

*M*

.

*Z m i*

*i*

*r*

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

## **2.3. Comments on multifractal and joint multifractal analysis**

## **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.

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

**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).

#### **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 � (scaling indices) and ���� (multifractal spectra), i.e., ���� and �������, corresponded to *q* = 15, and the maximum values, ���� and ������� corresponded to *q* = -15. A high value of ���� � ���� 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 ���� � ���� 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.

**3.1. Multifractal analysis**

the mapping of soil properties.

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

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 � (scaling indices) and ���� (multifractal spectra), i.e., ���� and �������, corresponded to *q* = 15, and the maximum values, ���� and ������� corresponded to *q* = -15. A high value of ���� � ���� 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 ���� � ���� 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

applications of multifractal and joint multifractal analysis in soil science.

**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 considered sufficient 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 (Ks) and organic carbon (OC) from the same study. For example, Fig. 4 shows the ���� (mass exponents) and ���� (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 Ks indicated multifractal behavior. Similarly, a large value for ���� � ���� of CL, OC, and Ks, 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].

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

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

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 ߙ 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 ߙ) 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

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

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

reported by Liu and Molz [46].

both water retention parameters, *n* and ߙ.

helped developing downscaling models [76].

**Figure 4.** a) Mass exponents of different variables. The solid line is the UM model passing through ��0�. 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].

122 Fractal Analysis and Chaos in Geosciences

reported with an increase in the EC value itself [63].

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 considered sufficient 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 (Ks) and organic carbon (OC) from the same study. For example, Fig. 4 shows the ���� (mass exponents) and ���� (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 Ks indicated multifractal behavior. Similarly, a large value for ���� � ���� of CL, OC, and Ks, 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

**Figure 4.** a) Mass exponents of different variables. The solid line is the UM model passing through ��0�. 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 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 ߙ 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 ߙ) 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, *n* and ߙ.

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.

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

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

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,

**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).

an abrupt textural change.

while the production of corn was monofractal [115].

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 (137Cs) fallout at small spatial scales has been found to be multifractal in nature [111]. Grubish [111] reviewed research on the 137Cs 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.

124 Fractal Analysis and Chaos in Geosciences

media [103].

the complexity of soil development at multiple scales.

characterize soil physical health and its quality [90].

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

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

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

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

The radioactive cesium (137Cs) fallout at small spatial scales has been found to be multifractal in nature [111]. Grubish [111] reviewed research on the 137Cs 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

simulated media for upscaling effective saturated hydraulic conductivity.

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).

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