search for


Multifractal Classification of the Disturbed Areas of the Sidi Chennane Phosphate Deposit, Morocco
Econ. Environ. Geol. 2022 Jun;55(3):231-9
Published online June 30, 2022;
Copyright © 2022 The Korean Society of Economic and Environmental Geology.

Abderrahim Ayad*, Saad Bakkali

Earth Sciences Department, Faculty of Sciences and Techniques, Abdelmalek Essaadi University, Tangier, Morocco
Received April 29, 2022; Revised June 7, 2022; Accepted June 7, 2022.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License ( which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
The irregular shape of the disturbances is a fundamental issue for mining engineers at the Sidi Chennane phosphate deposit in Morocco. A precise classification of disturbed areas is therefore necessary to understand their part in the overall volume of phosphate. In this paper, we investigate the theoretical and practical aspects of studying and measuring multifractal spectrums as a defining and representative parameter for distinguishing between the phosphate deposit of a low rate of disturbances and the deposit of a high rate. An empirical multifractal approach was used by analyzing the disturbed areas through the geoelectric images of an area located in the Sidi Chennane phosphate deposit. The Generalized fractal dimension, D(q), the Singularities of strength, α(q), the local dimension, f(α) and their conjugate parameter the mass exponent, τ(q) as well as f(α)-α spectrum were the common multifractal parameters used. The results reported show wide variations of the analyzed images, indicating that the multifractal analysis is an indicator for evaluate and characterize the disturbed areas within the phosphates deposits through the studied geoelectric images. This could be the starting point for future work aimed at improving phosphate exploration planning.
Keywords : phosphate series, local dimension, mass exponent, singularities of strength, generalized dimension
Research Highlights
  • The disturbances are one of the main issues encountered in Sidi Chennane phosphate deposit.

  • Geophysical mapping is used to detect the disturbed zones hidden under the quaternary cover.

  • The multifractal analysis is used to estimate the disturbances rate through the geophysical maps.

1. Introduction

Morocco is the world’s leading exporter of sedimentary phosphate with more than 30 percent of world exports and holds more than 50 billion MT of phosphate in reserve, which represents more than 70% of world supply. Phosphate rock has been known in Morocco since 1908, but its economic potential was not discovered until 1917 (Ayad and Bakkali, 2017; Brives, 1908). These phosphate reserves are sufficient to meet several centuries of demand. In this paper, the work was carried out in the Sidi Chennane deposit consisting of a succession of several high-quality phosphate layers. This deposit is part of the Ouled Abdoun basin and it is restricted to the west by the RP11 regional road, to the east by the RP309 national road, to the south by the 226000 meridian, and to the north by the outcrops of the lower phosphate layer (Fig. 1).

Figure 1. Geographical map of Sidi Chennane phosphate deposit located in the sedimentary basin of Ouled Abdoun. The study area is delimited by a red rectangle.

The phosphate series of Sidi Chennane frequently suffers from the inclusion of many sterile bodies called disturbances or derangements consisting of a waste rocky material and characterized with irregular geometrical shapes (Fig. 2)(Ayad et al., 2019; Ayad and Bakkali, 2019). Geological studies carried out in Sidi Chennane have revealed that the origin of these sterile bodies is linked to the dissolution of Senonian gypsiferous layers located under the phosphate series (Ait Taleb et al., 2009).

Figure 2. Section of phosphate sequences fronts demonstrating disturbances.

The presence of disturbances in the phosphate mines causes two major problems: (1) from an economic point of view, it is not easy to determine exactly the part of these bodies in the overall volume of the deposit and one cannot therefore correctly calculate phosphate reserves; (2) from a technical point of view, when encountering a disturbance, it is necessary to increase the number of holes in order to draw up the boundaries of the disturbances, and to increase the explosive power during the blasting operation, thereby explosives and time-consuming increase drastically.

Geophysical methods constitutes an appropriate tool in the Sidi Chennane deposit to cover the areas that may be disturbed (Bakkali, 2005). The anomalies detected can be used as a potential characteristic to create a model of the disturbances and would allow defining these structures before the mining front reaches them. In this work, we aim to numerically describe the visual aspect of the geoelectric images to diagnose them and quantify disturbed areas so that a control strategy can be developed for further estimation of phosphate reserve.

During these decades, many pattern recognition tools have been widely used to quantify complex structures. One of them is fractal and multifractal analysis used to describe the irregular geometry of many natural phenomena (Mandelbrot, 1989; Mandelbrot, 1990). In this paper, we intend to numerically characterize the shape of the disturbances using the following multifractal parameters: The generalized fractal dimension D(q), the singularities of strength α(q), the local dimension f(α), and their conjugate parameter the mass exponent τ(q) as well as f(α)-α spectrum. The analysis was performed on 2D geoelectric images of an area of 50 hectares located in the Sidi Chennane phosphate deposit.

2. Materials and Methods

2.1. Fractal and Multifractal

Fractal theory was invented by French mathematician Benoit Mandelbrot in 1975 (Mandelbrot, 1967; Mandelbrot, 1975). Mandelbrot introduced the fractal as an extremely irregular shape for which any suitably chosen part is similar to the whole at increasingly smaller scales, for which the Hausdorff-Besicovitch dimension exceeds its topological dimension as described in Mandelbrot et al. (1983).

The fractal approach has been widely used as an effective spatial analysis in different fields to better understand the fundamental scaling of many structures (Thomas et al., 2012; Tunas et al., 2016; Abidin et al., 2017, Ayad, 2022). Fractal analysis uses various techniques to assess and characterize the complexity of many objects using different fractal parameters such as fractal dimension (Ayad and Bakkali, 2022). However, the patterns of many complex objects cannot be fully described using simple fractal (monofractal) analysis. This analysis suffers from a lack of self-comparability and cannot be used to reveal much detail.

The weakness of fractal analysis can be corrected using multifractal analysis. It is a systematic approach that can help discover the characteristics of all kinds of objects and structures using more than one scaling process (Stanley et al., 1988). It provides more information about the space-filling characteristics than a fractal analysis. Compared to monofractal, a multifractal structure can be characterized as a superposition of homogeneous monofractal structures and cannot be described using the fractal dimension. Multifractal analysis has been widely employed successfully to characterize the complexity of many natural phenomena in different fields of scientific research (Lindinger et al., 2017; Mandal et al., 2017; Wang et al., 2017).

2.2. Geoelectrical Survey

Geoelectric prospecting is a tool for measuring the electrical resistivity of rocks below the surface. The working principle of this geophysical prospecting method is to inject an electric current of intensity I between two electrodes, and the voltage difference value (∆V) is measured between two potential electrodes after the current through the rocks. Then, the apparent resistivity of the rocks is found using the following formula: ρa = K (∆V/I), where K is the geometric factor of the instrument which depends on the configuration of the current and potential electrodes. The geoelectrical data used in this study was collected over an area of 50 hectares in the Sidi Chennane phosphate deposit. The measurements were performed by means of a Syscal2 resistivity meter using a 20m × 5m rectangular array. A total of 51 geoelectrical profiles spaced 20m apart were executed to reach a depth of 40m. There were 101 stations at 5m distance for each profile, which makes 5151 measurement stations.

The acquired data was then interpolated and imaged into a 2D contour map (Fig. 3). The map obtained allows defining the anomalies corresponding to the disturbances. It is used as a marker to distinguish between disturbances and phosphate rocks. The anomalies were classified as: (1) normal phosphate rocks of resistivity less than 200Ω.m; (2) disturbances of resistivity that mainly exceeded 200Ω.m.

Figure 3. Geoelectrical contour map of the study area.

2.3. Multifractal Analysis

To analyze the spatial pattern of disturbances, two sets of parameters are needed, including global and local parameters. The global parameters include the generalized correlation dimension D(q) and the mass exponent τ(q), while the local parameters include the singularity exponent (Lipschitz–Hölder exponent) α(q), and the local dimension (or Hausdorff dimension) f(α) of the set supporting this exponent (Chhabra and Jensen, 1989; Chhabra et al., 1989). The generalized correlation dimension D(q) is the basic multifractal parameter, based on the Renyi entropy. D(q) basically deals with the variation of mass varies as a function of ε (resolution or box size) in a digitized map. To calculate it, the spatial pattern of a fractal structure will be amplified step by step to some arbitrary exponent, q (-∞Hentschel and Procaccia (1983):


Where q is a continuous variable called the order of moment and Iq(ε) defines the Renyi entropy.

Iqε=1ni=1N εPi ε q

Where Pi(ε) is the distribution, indicating the ratio of the number of elements appearing in the ith box (rectangular grid) of linear scale ε to the number of the set, N(ε) represents the number of nonempty boxes.

The most frequently used generalized dimensions are D0 for q= 0; D1 for q= 1 and D2 for q= 2, which are respectively called the capacity dimension (or box counting dimension), information dimension (or entropy dimension) and correlation dimension.

The generalized dimension, D(q), is related to the other sets of multifractal exponents. Therefore, D(q), is obtained from the relationship with the second global parameter. τ(q) is known as the mass exponent function of q, which generalizes the relation given by:


The Hausdorff fractal dimension f(α) is a scale exponent of boxes with a common α, called Holder or singularity exponent. A graph of singularity spectrum can be plotted using f(α) versus α. The functions α and f(α) can be derived by the Legendre transformation as (Halsey et al., 1986):



The Legendre transform is mainly used to relate local parameters with global parameters. This suggests that the two sets of parameters can be converted into each other. If we calculate global parameters, we will obtain local parameters and vice versa.

In practice, from the τ(q) or D(q) curves, one can decide whether multifractality exists or not. If the curve τ(q) is a straight line or D(q) is a constant, then the structure is monofractal. However, if τ(q) or D(q) have convex shapes, then the structure is multifractal.

The multifractal process was performed using open source software called FracLac, a plugin for ImageJ (Karperien, 2004). The program uses the box-counting method to gather information about the distribution of pixel values in 2D binary images. Note that the box counting method consists to cover the image with a grid of decreasing box size, and then the boxes that include part of the figure are counted (Russel and Hanson, 1980). The collected data thus became the basis for generating graphs describing the multifractal scaling in order to show how the disturbances pattern behaves from one image to another.

In order to create a comparative analysis framework to obtain an accurate classification of the disturbed areas by the multifractal approach, we left sedimentary phosphates rocks out of our interpretation and we retain only the disturbed areas. Thus, the edges of the disturbances were traced and delineated by extracting their boundaries at different apparent resistivity cut-off frequencies (ʋar). Eight geoelectric images from ʋar=200Ω.m to ʋar=340Ω.m were created from the original image using a ʋa.r cut-off frequency of 20Ω.m. This allows observing how much the surface of the disturbances changes from one image to another. These maps were then converted to white and black binary images using Surfer 12 software. Where black pixels correspond to disturbed areas while white pixels represent sedimentary phosphates rocks (Fig. 4).

Figure 4. Spatial patterns of the disturbances extracted at eight cutoff frequencies through the geoelectrical contour map and their corresponding black and white images.
3. Results and Discussions

Before starting the analysis process, the main points of the multifractal approach must be explained. Multifractal spectrums have the potential to find differences in spatial information. It can provide us with multiple viewing angles from which to examine the surface more carefully. The global and local parameters were calculated for the order of moment -10 < q < 10. It can be noted that q can be used as a spatial analysis adjuster, by changing the value of q; we can analyze coarse structures and fine elements at different levels of view. The differences and similarities can be reflected by the multifractal spectra of the local and global parameters D(q), τ(q), α(q), and f(α), as shown in Fig. 5, 6, 7, and 8.

Figure 5. The generalized correlation dimension D(q) spectrums of the different images.
Figure 6. The mass exponent τ(q) spectrums of the different images.
Figure 7. The singularity exponent α(q) spectrums of the different images.

By using the global parameters (D(q) and τ(q)) (Fig. 5 and 6), we can analyze the whole spatial pattern of the disturbances. The D(q) curves are decreasing monotonic functions: D(q) ˃ D(q + 1). Roughly speaking, the functions D(q) are decreasing for negative values of q, while for positive values of q the curves are increasing. If q > 1, the high density pattern will be amplified, and therefore, the value of D(q) will decrease. Otherwise, if q < 1, the low density pattern will be amplified, and thus the value of D(q) will increase. In general, small differences in the disturbed surface may lead to significant changes in the resulting curves.

Moreover, for more precision, the simplest method for judging the results consists of comparing the capacity dimension, the information dimension, and the correlation dimension. Thus, since the spectrum D(q) is a monotonically decreasing function which decreases with increasing of the order of moment q, the capacity dimension D0 is necessarily greater than the information dimension D1 which is greater than the correlation dimension D2. Table 1 shows the results of D0, D1 and D2 of all analyzed images. It is obvious that, for all images, the inequality D0 > D1 > D2 is always established and respects the monotonically decreasing principle of D(q). The range of values (D0, D1 and D2) of the different analyzed images is between 0 and 2, i.e. 0

Table 1 . The calculated results of the capacity dimension, correlation dimension, and information dimension as a function of the corresponding rate of disturbances

ImagesDisturbances rate %Capacity dimension D0Information dimension D1Correlation dimension D2

The capacity dimension D0 is 1.558 for υ200Ω.m, and it decreases to 1.161 for υ340Ω.m, the information dimension D1 is 1.457 for υ200Ω.m, and it decreases to 0.995 for υ340Ω.m, and the correlation dimension D2 is 1.410 for υ200Ω.m, and it decreases to 0.916 for υ340Ω.m.

For the second global parameter, the mass exponent τ(q), show nonlinear curves. It is a monotonically increasing function τ(q)< τ(q + 1). Both left tails and the right ones of the curves show different slopes, indicating significant multifractal variability. The difference between these curves reveals the dissimilarity of the disturbed areas within the different images. In contrast, τ(q) tends to be linear for such pattern with lowest variation. In addition, when the statistical moment q = 1 the mass exponent is zero. In this case, there is a significant overlap of the curves. Moreover, it was noted that although the changes were not robust, these spectra could reveal significant differences.

By using the local parameters (f(α) and α(q)), we can analyze the small details of the disturbances (Fig. 7 and 8). This is different from the global parameters (D(q), τ(q)), which analyze the whole pattern. The singularity exponent curves α(q) show a decreasing monotonous behaviors α(q)˃ α(q + 1), while the spectra f(α) have the maximum f[α(q)]max=D0 when q = 0. Indeed, if q≤ 0, f(α) is a monotonically increasing quantity, and if q≥0, f(α) is a monotonically decreasing quantity. For q≤ 0, we have f [α(q)]≤f [[α(q + 1)] , while for q≥ 0, f [α(q)] ≥f [[α(q+ 1)]. The spectral behaviors of these parameters indicate a varying degree of multifractality for the analyzed images. It seems appropriate for estimating the degree of disturbances, which can allow a better classification of phosphate reserves (undisturbed).

Figure 8. The local dimension f(α) spectrums of the different images.

If the relationship between the local parameters f(α) and α(q) is plotted, an additional unimodal multifractal spectrum f(α)-α can be created (Fig. 9).

Figure 9. The singularity spectrums f(α)-α of the different images.

The f(α)-α spectra show an asymmetric shape, where the left tails are mostly shorter than the right ones. In all cases, these curves are found to be continuous with convex parabolas. When we comparing the spectrums f(α), we notice that their shapes are different, since all of the analyzed images are subject to different disturbances patterns. The dissimilarity of these curves indicates the possibility of associating these multifractal parameters to the rate of disturbances within the phosphate deposits.

Since the f(α)-α spectra does not provide a unique value, but a set of values, it is important to find the index which could be a discriminator in our application. For this purpose, the next step is to extract the indexes that describe these spectra. Several useful indices have been estimated in this work: αmin and αmax are respectively the minimum and maximum values of α, and f(αmin) f(αmax) (Table 2).

Table 2 . Results of the parameters extracted from the singularity spectrum as a function of the corresponding rate of disturbances

ImagesDisturbances rate %αminαmaxf(αmin)f(αmax)

According to Table 1, values of the estimated indexes can symbolize a discrimination of the analyzed images. The image ʋ200 which corresponds to the highest disturbance rate (17.683%) presents the highest values (αmin, αmin, f(αmin) and f(αmax)). However, the image ʋ340 which corresponds to the lowest disturbance rate (4.104%) presents the lowest multifractal values (αmin, αmin, f(αmin) and f(αmax)). We conclude that the increase in the surface occupied by the disturbances implies an increase in the variability of α, and vice versa.

4. Conclusion

In this research, we systematically analyze the potential utility of the multifractal approach to characterize the disturbed areas of the phosphate series of Sidi Chennane. Two sets of multifractal parameters are used: One is the global parameters, contain the generalized correlation dimension D(q) and the mass exponent τ(q), and the second is the local parameters include the singularity exponent α(q), and the Hausdorff fractal dimension f(α). The work was carried out on the geoelectric images of an area located in the northern part of the Sidi Chennane phosphate deposit, Ouled Abdoun, Morocco.

Eight geoelectric images were analyzed to make a comparative analysis and validate the potentiality of the multifractal approach. Clear results are obtained and suggest the strong sensitivity to characterize the spatial patterns of the disturbances through the geoelectric images. Comparing the similarities and differences can help determine how the surface of the disturbed areas changes, which may differentiate the phosphate deposit at high-risk of disturbances and the deposit at low-risk. The multifractal analysis approach can be used as an important analytical tool to make the best exploration and exploitation planning and obtain an automatic and accurate classification of phosphate reserves.

  1. Abidin, Ç. and Ulus, Ç. (2017) Three-dimensional modeling in medical image processing by using fractal geometry. J. Comput., v.12, p.479-485. doi: 10.17706/jcp.12.5.479-485
  2. Ait Taleb, Z., Mouflih, M., Benbouziane, A. and Amaghzaz, M. (2009) Description, pétrographie et origine des paléokasts du gisement de Sidi Chennane (Bassin des Oulad Abdoun, Maroc). Notes Mém. Serv. Géol. Maroc., v.530, p.21-30.
  3. Ayad, A. and Bakkali, S. (2022) Using the mass-radius method to quantify the disturbed zones in Sidi Chennane mine through geoelectrical images. International Journal of Mining and Geo-Engineering. doi: 10.22059/ijmge.2022.320175.594898
  4. Ayad, A. (2022) Prospecting for polymetallic mineralization in Middle Morocco using fractal mapping. Arab. J. Geosci., v.15, p.1-12. doi: 10.1007/s12517-022-09613-2
  5. Ayad, A. and Bakkali S. (2019) fractal assessment of the disturbances of phosphate series using lacunarity and succolarity analysis on geoelectrical images (Sidi Chennane, Morocco). Complex., v.2019, p.1-12. doi: 10.1155/2019/9404567
  6. Ayad, A., Amrani, M. and Bakkali, S. (2019) quantification of the disturbances of phosphate series using the box-counting method on geoelectrical images (Sidi Chennane, Morocco). Int. J. Geophys., v.2019, p.1-12. doi: 10.1155/2019/2565430
  7. Ayad, A. and Bakkali, S. (2017) interpretation of potential gravity anomalies of Ouled Abdoun phosphate basin (Central Morocco). J. Mater. Environ. Sci., v.8, p.3391-3397.
  8. Bakkali, S. (2005) Analysis of phosphate deposit “disturbances” using the horizontal-gradient responses of resistivity data (Oulad Abdoun, Morocco). Earth Sci. Res. J., v.9, p.123-131.
  9. Brives, A. (1908) Sur le Sénonien et l’Eocène de la bordure nord de l’Atlas marocain. C. R. Acad. Sci. Ser. Paris., v.196, p.873-875.
  10. Chhabra, A. and Jensen, R.V. (1989) Direct determination of the f(α) singularity spectrum. Phys. Rev. Lett., v.62, p.1327-1330. doi: 10.1103/physrevlett.62.1327
    Pubmed CrossRef
  11. Chhabra, A., Meneveau, C., Jensen, R.V. and Sreenivasan, K.R. (1989) Direct determination of the f(α) singularity spectrum and its application to fully developed turbulence. Phys. Rev. A, v.40, p.5284-5294. doi: 10.1103/physreva.40.5284
    Pubmed CrossRef
  12. Halsey, T.C., Jensen, M.H., Kadanoff, L.P., Procaccia, I. and Shraiman, B.I. (1986) Fractal measures and their singularities: The characterization of strange sets. Phys. Rev. A., v.33, p.1141-1151. doi: 10.1103/physreva.33.1141
    Pubmed CrossRef
  13. Hentschel, H.G.E. and Procaccia, I. (1983) The infinite number of generalized dimensions of fractals and strange attractors. J. Phys. D., v.8, p.435-444. doi: 10.1016/0167-2789(83)90235-x
  14. Karperien, A. (2004) FracLac Advanced User’s Manual (Frac_Lac.jar) Plugins/Fractal Analysis/FracLac. 2005/03/23. version 2.0aF. Charles Stuart University, Australia. Available online at
  15. Lindinger, J. and Rodríguez, A. (2017) Multifractal finite-size scaling at the Anderson transition in the unitary symmetry class. Phys. Rev. B., v.96, p.1-11. doi: 10.1103/physrevb.96.134202
  16. Mandal, S., Roychowdhury, T., Chirom, K., Bhattacharya, A. and Brojen, S.R.K. (2017) Complex multifractal nature in Mycobacterium tuberculosis genome. Sci. Rep., v.7, p.1-13. doi: 10.1038/srep46395 Mandelbrot, B.B. Freeman, W.H. and San Francisco, C.O. (1983)
    Pubmed KoreaMed CrossRef
  17. The Fractal Geometry of Nature. Earth Surf. Process. Landf., v.8, p.460. doi: 10.1002/esp.3290080415
  18. Mandelbrot, B.B. (1967) How Long Is the Coast of Britain? Statistical Self-Similarity and Fractal Dimension. Sci., v.156, p.636-638. doi: 10.1126/science.156.3775.636
    Pubmed CrossRef
  19. Mandelbrot, B.B. (1989) Multifractal measures, especially for the geophysicist. Pure Appl. Geophys., v.131, p.5-42. doi: 10.1007/bf00874478
  20. Mandelbrot, B.B. (1990) Negative fractal dimensions and multifractals. Phys. A., v.163, p.306-315. doi: 10.1016/0378-4371(90)90339-t
  21. Mandelbrot, B.B. (1975) Stochastic Models of the Earth's Relief, the Shape and the Fractal Dimension of the Coastlines, and the Number-area Rule for Islands. Proc. Natl. Acad. Sci. U.S.A., v.72, p.3825-3828. doi: 10.1073/pnas.72.10.3825
    Pubmed KoreaMed CrossRef
  22. Russel, D., Hanson, J. and Ott, E. (1980) Dimension of strange attractors. Phys. Rev. Lett., v.45, p.1175-1178. doi: 10.1103/physrevlett.45.1175
  23. Stanley, H.E. and Meakin, P. (1988) Multifractal phenomena in physics and chemistry. Nature, v.335, p.405-409. doi: 10.1038/335405a0
  24. Thomas, I., Frankhauser, P. and Badariotti, D. (2012) Comparing the fractality of European urban neighbourhoods: do national contexts matter. J. Geogr. Syst., v.14, p.189-208. doi: 10.1007/s10109-010-0142-4
  25. Tunas, I.G., Anwar, N. and Lasmint, U. (2016) Fractal Characteristic Analysis of Watershed as Variable of Synthetic Unit Hydrograph Model. Open Civ. Eng. J., v.10, p.706-718. doi: 10.2174/1874149501610010706
  26. Wang, W., Cheng, Q., Tang, J., Pubuciren, Song, Y., Li, Y. and Liu, Z. (2017) Fractal/multifractal analysis in support of mineral exploration in the Duolong mineral district, Tibet, China. Geochem.: Explor. Environ. Anal., v.17, p.261-276. doi: 10.1144/geochem2016-449


June 2022, 55 (3)