The phosphate series of Sidi Chennane (Morocco) is frequently affected by sterile bodies called “disturbances” causing many problems during the mining operation. The quantification of these bodies is thus considered as a crucial step for the OCP mining engineers. In following this research, we propose the fractal dimension as a new and simple efficient tool to analyze the disturbed areas. The work was carried out on geoelectrical maps of a study area of 50 hectares located in the northern part of the Sidi Chennane deposit. Fractal dimension was used as a representative parameter for examining the disturbed areas by using the box-counting method. The end result shows a strong linkage between the rate of the disturbances and the corresponding fractal dimension. This may have an important implication to quantify the disturbances and get a precise phosphate reserves estimate.

1. Introduction

Morocco holds more than 75% of the global phosphate reserves with more than 35 billion cubic meters. It is the world’s leading exporter and the third largest producer after the USA and China with an annual output of 19 million tons. These appreciable reserves of phosphate were deposited in several basins, from northeast to southwest of Morocco (Ouled Abdoun, Gantour, Meskala, and Boukraa) and have been explored by L’Office Chérifien des Phosphates (OCP) since the 1920s [1].

In this study, we focused on Sidi Chennane mine which is one of the most important phosphate deposits in Ouled Abdoun sedimentary basin [2], located near the city of Khouribga, and at about 100 km in the southeast of Casablanca (Figure 1).

The phosphate series of Sidi Chennane frequently contains many sterile bodies of extremely locally high density called disturbances as a waste rocky material that must be removed (Figure 2).

These bodies are characterized by an irregular geometrical shape and dispersed without any order or privileged direction under a Quaternary cover. The geological investigations carried out in Sidi Chennane deposit have revealed that the origin of these structures is linked to the dissolution of the Senonian gypsiferous layers situated at the base of the phosphate series [3]. Their detection and delimitation required in many cases several exploration works, such as drilling wells that remain not particularly effective and increase drastically the cost of phosphate extraction. In light of this issue, the OCP Group has implemented with some Moroccan universities the geophysical mapping methods that underline the contrast of density between the disturbances and normal phosphate layers [4, 5]. This had an important implication to delineate the disturbed areas in order to get a precise phosphate reserves estimate.

In the exploration and the estimation of the phosphate reserve, we need to know the local surface of the disturbances. The OCP Group engineers usually use classical approaches including geostatistical and geometrical methods. The disturbed areas are partitioned according to the exploration plan and the phosphate reserves will be thus estimated through a discrete and linear function. Using these classical geometric methods, the irregular surface of the disturbances is only roughly estimated. The same surface measured at different scales will give different results.

In the last two decades, fractal geometry gained large popularity essentially in geosciences, where the measures to describe the forms of surfaces are essential. It was applied extensively in many studies and applications like rocks texture analysis [6], natural fractures pattern recognition [7], porous media [811], and naturally fractured reservoirs [1214]. Throughout the present work, we have tried to point out that fractal geometry seems promising and more appropriate to quantify the disturbed areas since they commonly show an irregular geometry on the geoelectrical maps [15].

To achieve the goal of this paper, we present and discuss a case study carried out on an area of 50 hectares located in Sidi Chennane deposit. The applicability and the relevance of fractal geometry in this study are verified from the variability of the geometrical aspect of the disturbances on the geoelectrical maps. This study aims to provide quantitative and useful means to quantify the disturbances and understand how going back fractal geometry may help in the phosphate reserve estimation.

1.1. Fractal Geometry Background

The fractals were coined by Benoît Mandelbrot in 1974 [16]. The term “fractal” starts from the Latin root “fractus” perfect passive participle of “frang” which means an irregular structure created following stochastic or deterministic rules. The fractal geometry is used to describe objects of complex shape that are practically impossible to be evaluated by conventional methods based on Euclidean geometry. Mandelbrot began his treatise on fractal geometry by explaining that a wide variety of natural objects does not conform to Euclidean geometry, “a mountain is not a cone, a coastline is not a circle, and clouds are not spheres”[17]. For example, if we ask “How long is the coast of Britain?,” the answer to this question depends on how the measuring stick length is. As we let the length of the stick get small, the length of the coast becomes infinite, because the surface of the coast shows increasing detail the closer we zoom into it. A fractal is thus an object which consists of distinctive elements generated by the process of repeating the same operation for retaining the same pattern at different scales of examination. This is the central concept of fractal geometry summarized in the notion “self-similarity.”

Over the years, fractal analysis becomes well known as a valuable tool for a quantitative measure of surface roughness. It is introduced as being complementary to solve the complex character of objects whose shapes show an irregularity. The main attraction of this approach stems from its ability to offer an appropriate way to quantitatively analyze data in a variety of fields, such as economics, physics, earth sciences, biology, chemistry, and other disciplines [19, 20]. In contrast to Euclidian geometry measurement, fractal geometry is a relative way that allows obtaining the dimension of a very rough or bumpy area, length, or volume. In the Euclidian geometry, each set of elementary geometry has an integer dimension that equals 1 for straight line, 2 for square or circle, and 3 for the cube. This intuitive dimension is commonly known as Euclidean or topological dimension. However, fractal geometry allows the complex or irregular objects to be characterized in noninteger dimensions with values that may vary anywhere between 1 and 2 called fractal dimension FD. Figure 3 illustrates an excellent representation of the fractal dimension variation from 1 to 2D.

Moreover, many mathematical functions are considered space filling which maps a one-dimensional interval into a two-dimensional area such as Cantor set, Sierpiński curve, Dragon curve, and Peano curve [2124]. These fractals are called mathematical fractals generated by exact iteration and do not have scaling limits.

2. Materials and Methods

Many methods for measurement fractal dimension can be used and the results obtained by different methods often differ significantly. Furthermore, not only the method but also the software used to calculate the fractal dimension may contribute to the differences. In this study, we obtained the fractal dimension by using the most popular method called the box-counting method. It is the well-known fractal method that is proposed by Russel [25]. It is an easier tool in the calculation process and more effective to deal with the irregularity. The typical approach of this method consists in covering the image to be measured with a grid of decreasing box size r and then counting the boxes elements that include some portion of the object N(r) (Figure 4).

The goal is to evaluate how the size of the box used to inspect pixels affects the number of boxes it takes to inspect all of the pixels. From one step to the next, as the size of the boxes shrink, the number of occupied boxes increases. Therefore, we obtain a set of data points that will be represented on a log-log Cartesian graph, which allows for identifying the fractal dimension value through a biasymptotic curve (named the empirical curve). For typical fractal object examined using the box-counting approach, the probability of a number of boxes N(r) includes some portion of the figure that is defined following (1), with FD being the fractal dimension:The number of boxes (N) and box size (r) are expressed as a power law. The fractal dimension will be thus defined as where r is the constant size of the box and N(r) is the number of boxes needed to completely cover the object in each measurement. The image is thus considered fractal when the relationship between the logarithms of N(r) versus the logarithms of (r) is linear on a graph. This linearity indicates that a power law relates the number of boxes covering the object to their size.

In this work, the fractal analysis was performed on a set of 5000 apparent resistivity data acquired from a study area of 50 hectares within Sidi Chennane deposit [26, 27]. The use of the geoelectrical method has its main aim in the detection of potential anomalies, and it is used in this work as an excellent parameter and a marker for distinguishing between disturbances and the phosphate rocks [28]. The geoelectrical anomalies are classified into phosphate rocks of a resistivity less than 200Ω.m, and the disturbances that correspond to high values of apparent resistivity exceeded 200Ω.m (Figure 5).

To get a complete characterization of the disturbed areas and provide accurate information, a set of fractal dimensions is needed, showing how the different disturbed areas have different fractal dimensions. For this aim, we start by constructing eight maps of apparent resistivity range from 200Ω.m to 340Ω.m using a cutoff frequency ʋa.r of 20Ω.m (Figure 6).

The box-counting analysis was carried on in this study by means of two separate programs HarFA 5.5 and fractalyse 2.4. Considering that using the combination of these two software programs would serve to cross-check and validate the results. These software programs were both used by researchers worldwide to approximate the complexity of objects, and they seem to be reliable as they have been validated in many studies [29, 30].

Before starting the analysis process, the geoelectrical maps have to be transformed into black and white binary images where the black pixels correspond to disturbances and white pixels correspond to phosphate rocks (Figure 8). The images were then transformed on bitmap format in order to be adapted to the fractal analysis software requirement.

Quantitative observations about geoelectrical images can be made to show a view of the geometric forms of disturbed areas. The edges extracted at different cutoff frequencies ʋa.r allow observing to what extent the surface of the disturbance changes from one image to another. Then, for each image, the surface of the disturbances was explored by means of fractal dimension using the box-counting method.

Moreover, the precise estimation of the fractal dimension requires the choice of adequate size for the analyzed images. Because the image size is one of the most important parameters which impacts the results, changing the size may give different results for the same image. In light of this issue, the fractal dimension was firstly calculated for an image with full filling space (black image) without any irregularity or roughness, which corresponds generally to 2D. This is considered as an experimental basis for the irregular character of the disturbances. In order to obtain the correct result, we carried out several tests using different sizes and we have obtained different values of the fractal dimension. The chosen size corresponds to 510 pixels by 255 pixels which gave the best fit of fractal dimension, the closest to 2 (1.9995 by means of HarFA 5.5 with an error of ±0.0002 and 1.998 by means of fractalyse 2.4 with an error of ±0.0008) (Figure 7). Furthermore, as one can see the fractal analysis programs (HarFA 5.5 and fractalyse 2.4) gave slightly different results, which suggested that these software programs were operating correctly and were accurately characterized the fractality of the image submitted.

From this result, it is obvious that the fractal dimension FD≈2 would correspond to an area without any disturbances. Meanwhile, this value will decrease and will be in the range of 1-2 according to the percentage of disturbances in which each one of the images is illustrated in Figure 7.

3. Results and Discussion

The validity of this study has been thus verified by applying the fractal analysis to images of different disturbed surfaces. Then, the number of boxes computed N(r) were plotted in a typical ln-ln graph (double logarithmic plots) versus their side length (r). The graphs are depicted in both Figures 9 and 10.

It is shown in Figures 9 and 10 that the fractal dimension computed exhibits a dependence between the logarithm of the number of boxes necessary to cover the disturbances and the logarithm of the box sizes. These plots prove different points, the most important of which is that the surface of the disturbances proves, in fact, to be fractal in a mathematical point of view. All the maps analyzed exhibited strong fractal properties, as one would note the linear character of the relationship between the logarithms of the number of occupied boxes and the box sizes highlighted by the high values of the correlation coefficients R2. Moreover, it is important to note that employing both two software programs in measurement added useful detail to cross-check and validate the results.

In this spatial analysis, it is interesting not only to characterize the morphology of the disturbances but also to help get an overall quantitative understanding of how the disturbed surface varies as a function of the fractal dimension. For this aim, it is important to systematically classify and compare between them. Firstly, the surface of the disturbances was measured in each map using a Geographic Information System (GIS). Afterwards, the fractal dimensions calculated from both programs were significantly correlated with their rates of disturbances for showing a clear view of the analysis. The results obtained are summarized in Table 1.

From Table 1, it seems that the fractal dimension grows constantly with the increase of the surface of disturbances. The fractal dimensions calculated are in the range of 1.40-1.66 by means of HarFA 5.5 () and in the range of 1.39-1.62 by means of fractalyse 2.4 (). For each one of the software programs, the FD values correspond, respectively, to the rates of disturbances range from 4.1% to 17.7%. While this must emphasize that the fractal dimension using fractalyse 2.4 showed smaller values than the fractal dimension using HarFA 5.5, this divergence is more pronounced in maps having the highest disturbances rate.

It is apparent that there is an important statistical correlation that allows an expeditious comparison and thus may help in decision making as an index in phosphate exploration planning. The variation of fractal dimension values seems reasonable from viewing the different sampled images (Figure 7). It appears that the variation of FD values was really influenced by the spatial distribution of the disturbances and their shape into the overall area. This observation provides proof to support the process that a fractal dimension is a measure of space-filling capacity. As the space-filling rate is increased the fractal dimension is also increased and vice versa, indexing how the detail in disturbed areas changes. It reflects the relative amount of variation in the pattern, such that patterns with higher space-filling and thus higher visual complexity are quantified by higher FD values. Hence, the highest values of FD indicated that the disturbed areas are filling more space, resulting thus in a more shattered pattern. In contrast, the lowest fractal dimensions (FD) are typically associated with disturbed areas filling less space. This positive correlation indicates that the fractal dimension is an effective measure that proves this significance and seems to be most appropriate to quantify, compare, and evaluate the rate of the disturbances when applied to geoelectrical images.

The present study has improved the way of phosphate exploration planning through a comprehensive understanding of the spatial character of the disturbances. This is often a major concern in the Sidi Chennane area and other deposits around Ouled Abdoun basin. This analysis has described how much space is filled by disturbances in the area wherein they are embedded and therefore can help in the phosphate reserves estimation. Furthermore, the study is simple enough to be applicable to any phosphate deposits in Ouled Abdoun, as such opening a window of opportunity for the development of other automatic and quantitative tools to estimate the rate of the disturbances.

4. Conclusions

Throughout this work, we studied the fractal dimension utility to quantify the disturbances of the phosphate series on the geoelectrical maps. The box-counting method used to calculate the fractal dimension, and its performances were described. The article brings also the principle and the methodological frame of this method including formulas by which the fractal dimension was calculated.

Detailed examination of relations between disturbed areas and their corresponding fractal dimensions was achieved by means of two different software programs, HarFA 5.5 and fractalyse 2.4, which provides an interesting comparative and contextual framework for analysis. The end results suggested confirm the goodness for using the fractal dimension to examine the spatial patterns of the disturbances.

Since the OCP Group managers are eager for developing tools to quantify the impact of disturbances, the fractal analysis has met management needs since it produces results of practical significance. This framework of analysis may lead to an interesting assumption about using fractal geometry to get accurate phosphate reserve estimation and make the best exploration planning.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this work.


The authors thank the anonymous reviewers of the International Journal of Geophysics. This study was performed at the Environment, Oceanology, and Natural Resources Laboratory, Faculty of Sciences and Techniques of Tangiers, Morocco. The research is supported by a doctoral grant [UAE, L002/005, 2016] from the National Center for Scientific and Technical Research of Morocco (CNRST).