Printer Friendly

Geostatistical analysis for delineating sterile inclusions in Sidi Chennane' phosphatic series, Morrocco.


Morocco is the world's third largest phosphate producer, after the USA and China. Total mine production recorded by the Ministry of Energy and Mines in 2013 was 30 Mt so more than 75% of world reserves. Four major phosphate basins are now known and are being exploited, three of which are located in central-northern Morocco. The four main deposits of phosphate are: the Oued Eddahab basin situated in Sahara, the central ganntour basin near Youssoufia, the Meskala basin at east of Essaouira and the Oulad Abdoun basin situated near Khouribga. The existence of morocco sedimentary phosphate rock has been known since 1908 in the Meskala basin, but it had not generated significant interest until the discovery, in 1917, of the Oulad Abdoun basin.

The geological investigations carried out in Sidi Chennane phosphatic deposit in the Oulad Abdoun basin revealed a phenomenon of sterile hardpan [Bakkali 2005, Baba & al. 2012; 2014a and 2014b] these sterile bodies are formed by accumulations of silicified limestones or by limestone blocks within an argillaceous matrix which interfere with phosphate extraction and their resistivity, above 200 Ohm.m, is higher than the phosphate-rich mineral resistivity wich is near 120 Ohm.m [Bakkali and Amrani 2008; Baba & al. 2012]. The application of the electric prospection methods constitutes a suitable means to map these sterile bodies [Baba & al. 2012, 2014a and 2014b].

One of tools for mapping of electrical resistivity data is geostatistical methods. Resistivity data were analyzed using geostatistics (variogram) in order to describe and quantify the spatial continuity.

This study has been carried out with the aim of evaluation geostatistical methods for estimation of subsurface layers delimitation and specified their spatial variability in order to establish a model of sterile bodies' distribution and would permit the definition of these structures before the mining front reaches them.


Data acquisition

Electrical and electromagnetic (EM) methods have been important in the field of Applied Geophysics for about a century, particularly for shallow and near-surface investigations. Geoelectrical resistivity surveying has become an important and useful tool in hydrogeological studies, mineral prospecting and mining, as well as in environmental and engineering applications [Griffiths & al. 1990; Griffiths and Barker, 1993; Dahlin and Loke, 1998; Olayinka, 1999; Olayinka and Yaramanci, 1999; Amidu and Olayinka, 2006; Baba & al. 2012; Ouadif & al. 2012].

Electrical sounding is a method to investigate the change in earth resistivity with depth at a particular location. Traditionally, arrangements using four electrodes (two current-transmitting electrodes and two voltagesensing electrodes) are used for either vertical soundings or horizontal profiling. For vertical soundings, the electrodes are arranged symmetrically according to a center, with increasing distances between electrodes used to explore deeper depths [Baba & al. 2012; Ouadif & al. 2012 and 2014a].

When using geoelectrical surveying techniques, an estimate of the earth resistivity is calculated by using the well-known relation between resistivity, an electric field, current density (called Ohm's Law), the geometry and spacing of the current and potential electrodes. When the earth is not homogeneous neither isotropic, this estimated value is called the apparent resistivity, p_app, which is an average of the true resistivity in the measured section of the earth.

Profiling resistivity surveys can be employed to detect lateral variations in resistivity. This type of survey provides estimates of the spatial variation in resistivity at some fixed electrode spacing [Chapellier 2000; Dabas and Caraire 2012].

To determine the intrinsic resistivity for the different terms of the phosphatic series, we carried out, during the geophysical prospection in a parcel of 50 ha,5151 resistivity measurements as horizontal profiling (SYSCAL-R2 resistivity instrument IRIS Instruments) using the wellknown Schlumberger array, in order to map the spatial distribution of the sterile hardpan inclusions. The 5151 stations of the resistivity survey is a compilation of 51 profiles spaced at 20 m. For every profile, there were 101 stations separated by 5 m using the well-known Schlumberger array.

1.1 Geostatistics

Geostatistics, started with the work in geology and mining of Krige (1951) then developed by Matheron (1963) with his theory of regionalized variables, are methods of interpolation which predict unknown values from data observed at known locations, and it minimizes the error of predicted values which are estimated by spatial distribution of the predicted values. As variables are supposed having an effect on each other, the first concept is related to the variogram of values we want to estimate. It relates the similarity or difference, expressed as the semi-variance, between values at different places to their separation distance (lag) and direction. Distance that variables have an effect on each other is called range of a variogram and is an important concept because after that distance, variables have no effect on each other.

The plot of the semivariance versus the lag h is called semivariogram. In some cases, it increases with h and reaches an asymptotic value, a sill. Range is the value of h at which the semivariogram asymptotes. Values of the properties at sites separated by a distance shorter than the range will be spatially correlated. The range depends on the scale of the study and on the property evaluated.

positions x and x + h, with h a distance, known as la g, which separates both positions, semi-variances can be estimated as:

[gamma](h) = 1/2 E [[[rho](x) - [rho][x + h]].sup.2] (1)

For discrete sampling sites, such as earthen material samples, the function is estimated as:

[gamma](h) = 1/2N(h) [N(h).summation over (i=1)] [[rho]([x.sub.i]) - [rho][[x.sub.i] + h]].sup.2] (2)

Where [rho]([x.sub.i]) : the value of electrical resistivity [rho] at location of [x.sub.i]

N(h) : the number of point pairs for distance class h.

The five models most frequently used are:



Exponential: g(h)=c. (1 - exp (-3h/a)) (4)

Gaussian: g(h)=c. (1 - exp (-3[h.sup.2]/[a.sup.2])) (5)

Power: g(h)=c.[h.sup.w] o<w<2 (6)

With h: lag distance

c: sill

a: range

Nugget effect is essential to optimize the process under study: the higher its difference for the semivariogram sill, the higher the phenomenon continuity, the lower the estimation variance, or the higher the estimation confidence.

The variogram plot is fitted with a theoretical model which provides information about the spatial structure as well as the input parameters for kriging interpolation. Kriging algorithms (Krige, 1951) are described as a best linear unbiased estimator (BLUE), which is a process of a theoretical weighted moving average:

[rho]([x.sub.0]) = [n.summation over (i=1)] [[lambda].sub.i][rho]([x.sub.i]) (7)

where [rho]([x.sub.0]) is the value to be estimated at the location of [x.sub.0], [rho]([x.sub.i]) is the known value at the sampling site xi. And [[lambda].sub.i] denotes the weight of the [] observation, they are calculated searching an unbiased estimator of z with minimum variance [Webster, 1985].

The estimation errors (or kriging variances) need to be minimised. To ensure that the estimate is unbiased, the weights need to sum to one:

[n.summation over (i=1)] [[lambda].sub.i] = 1 (8)


1.1. Geological framework and methodology of work

The Ouled Abdoun basin is the largest phosphate basin in Morocco. It's located about 100 km in south-east of Casablanca. The phosphate deposits of Ouled Abdoun area belongs to the western Moroccan Meseta, commonly considered being stable. The local sedimentary deposits resulting from a large transgression occurred in mid-Cretaceous. It consists of :marly limestone and gypsum of Cenomanian, Turonian white limestones, Senonian marl and yellow marly limestones, phosphatic series dated from Maastrichtian to Ypresian and Lutetian calcareous Thersitean slab. The Neogene continental deposits cover locally the marine series.

1.2 Statistical Analysis

The classical statistical procedures assume that variations are randomly distributed within each class of observations and independent of sample location [Mohmmadi, 2006]. The mean is therefore used for the estimation of properties for unsampled locations within each class. Statistics of dispersion (e.g. variance, standard deviation, coefficient of variability, confidence limits) are used to show the precision of the mean as an estimator.

The resistivity distributions are not normal, which is the particularity of the anomalous bearing population. Variation between maximum and minimum for these data shows a wide range.

Assuming the samples were collected randomly from different locations in the field, we can conclude from these statistics that resistivity at the upper sublayers is more variable on the field than that at the lower layers.

2.3 Geostatistical methods

Vesper 1.6 software (ACPA, University of Sydney) developed by the Australian Center for Precision Farming [Minasny & al. 2002] was used to calculate experimental electrical resistivity semivariograms of experimental fields in order to obtain the parameters: range, sill and nugget variance for evaluating the spatial correlation of these variables. Semivariograms were obtained using a classical Matheron semivariogram estimator.

The experimental semivariograms were fitted with various theoretical models like spherical, exponential, gaussian, linear and power by an automatic (least-squares) estimation of semivariogram parameters. Indicator variograms were built adopting a lag of 50 m. The theoretical model that gave minimum standard error is chosen for further analysis.

Table 2 shows semivariogram models and parameters. The attributes were fitted to the exponential model.

The ratio of nugget to total semivariance, expressed as a percentage, was used to classify spatial dependence [Cambardella & al. 1994], a ratio less than 25% indicated strong spatial dependence, between 25 and 75% indicated moderate spatial dependence and more than 75% indicated weak spatial dependence.

Spatial interpolation accuracy and precision were evaluated through crossvalidation approach. Most important criteria include: mean bias error (MBE), mean absolute error (MAE), mean squared error (MSE) and root mean squared error (RMSE). MAE and RMSE are among the "best" overall measures of model performance [Willmott, 1982, Webster and Oliver 2001]. RMSE provides a measure of error size, but it is sensitive to outliers as it places a lot of weight on large errors [Chai and Draxler 2014]. MAE provides an absolute measure of the size of the error and it is less sensitive to extreme values [Delbari & al. 2013]. The definitions of MAE, MBE and RMSE are as follows:

MAE = 1/n [[SIGMA].sup.n.sub.i=1][absolute value of [z.sup.*]([x.sub.i]) - z([x.sub.i])] (9)

MBE = 1/n [[SIGMA].sup.n.sub.i=1] ([z.sup.*]([x.sub.i]) - z([x.sub.i])) (10)

RMSE = [[square root of 1/n [[SIGMA].sup.n.sub.i=1] ([z.sup.*]([x.sub.i]) - z([x.sub.i]))].sup.2] (11)

Where z(xi) and [z.sup.*](xi) are the observed and predicted value respectively. To see how predicted values match closely the observed values, the RMSE can be used: the smaller the RMSE the better.

We used the RMSE to evaluate model performances in cross-validation mode. The smallest RMSE indicate the most accurate predictions.

Cross-validation of the model are shown in figures 6.a, 6.b and 6.c for the variables AB40, AB80 and AB120. Because the RMSE are closer to zero, fitted Exponential models yield better results and are adopted in this study.

After analysis of spatial dependence, we evaluated sites that had not been sampled in order to produce surface maps by ordinary kriging.

Figures 7.a, 7.b and 7.c show isoresistivity Maps carried out at the nodes of the square grid of 1km x 1km with the Exponential model of variogram.

The comparison of the maps of the phosphate deposit anomalous zones corresponding to the three Schlumberger devices improved the quality of subsurface layers delimitation and specified their spatial variability: (1) for AB=40 and AB=80 the contour maps show anomalous zones located at an average depth of 20m (figures 7.a and 7.b). (2) For AB=120 (figure 7.c) the disturbances as detected from surface measurements, which present depth higher than 20 m, seem to be not very significant regarding to surface study area. Thus, the majority of anomalies affect only the subsurface of the phosphate series


In this work ordinary kriging, a type of geostatistical techniques is applied for mapping apparent resistivity at Sidi Chennane, one of phosphatic basins in Morocco, which is characterized by the presence of sterile hardpan -so-called "disturbances"- hard to detect. A detailed geoelectric resistivity survey was executed using the horizontal electrical sounding approach with shlumberger array, in order to study the lateral variations in the geoelectric behavior.

The data process was executed using the Vesper 1.6 software (ACPA, University of Sydney). The Exponential model is found to be the best model representing the spatial variability of the resistivity, for this reason, it is used to plot resistivity map by kriging method. These maps permitted us to identify disturbed zones in this area; the qualitative interpretation of these maps reflects that the disturbances as detected from surface measurements which present depth higher than 20 m seem to be not very significant.

The study showed the benefits of the using geostatistical methods in the assessment of data, sampling strategy, and estimation of values at unsampled areas of the sub-surface to delineate anomalies. Our results suggest that geostatistical methods are potentially very useful for building a model in the area of complicated geological structures.


Amidu S.A., Olayinka A.I. 2006. Environmental assessment of sewage disposal systems using 2D electrical resistivity imaging and geochemical analysis: A case study from Ibadan, Southwestern Nigeria. Environ. Eng. Geosci., Vol.7, No. 3, pp. 261-272.

Baba, K., Bahi, L., Ouadif, L. & Akhssas, A. 2012. Mapping Sterile Bodies in the Sidi Chennane phosphatic deposit (Morocco) using geoelectrical Investigations. International Journal of Engineering Research and Applications (IJERA), Vol. 2, No. 5, pp. 2132-2136.

Baba, K., Bahi L. and Ouadif, L. 2014a. Enhancing geophysical signals through the use of Savitzky-Golay filtering method. Geofisica Internacional. Vol.53, No. 4: pp. 153-162.

Baba, K., Bahi, L., Ouadif, L. & Cherradi, C. 2014b. Application des methodes d'analyses statistiques multivariees a la delimitation des anomalies de Sidi Chennane. Journal des materiaux et des sciences de l'environnement Vol.5, No. 4, pp. 1005-1012.

Bakkali, S. 2005. Analysis of phosphate deposit disturbances using the horizontal gradient responses of resistivity data (Oulad Abdoune Morocco). Earth Sci. Res. J., Vol. 9, No. 2, pp. 123-131.

Bakkali, S. and Amrani, M. 2008. Denoising resistivity phosphate "disturbances" using haar mother wavelet transform (Sidi Chennane, Morocco). Earth Sci. Res. J., Vol. 12, No. 1, pp. 62-71.

Cambardella, C.A.; Moorman, T.B.; Novark, J.M.; Parkin, T.B.; Karlen, D.L.; Turco, R.F.; Konopka, A.E. (1994). Field-scale variability of soil properties in central Iowa soils. Soil Science Society of America Journal 58, pp. 1501-1511.

Chai, T. and Draxler, R. R. 2014. Root mean Square error or mean absolute error? Arguments against avoiding RMSE in the literature. Geoscientific Model Development, 7, pp. 1247-1250.

Chapellier, D. 2000. Prospection Electrique en Surface. Cours de Geophysique. Universite de Lausanne, Institut Francais de Petrole.

Dabas, M. and Caraire, G. 2012. Arp Imaging for vulnerability assessment of RTE Grounded high voltage power line. 100 Years of Electrical Imaging, Collection Sciences de la Terre et de l'environnement, Presses des Mines, pp. 29-33.

Dahlin, T. and Loke M.H. 1998. Resolution of 2D Wenner resistivity imaging as assessed by numerical modelling. J. Applied Geophys., 38, pp.237-249.

Delbari M., Bahraini M. M. and Amiri M. 2013. Spatio-temporal variability of groundwater depth in the Eghlid aquifer in southern Iran. Earth Sci. Res. J. vol.17 No. 2, pp. 105-114

Griffiths D.H., turnbull, J. and olayinka, A.I. 1990. Two-dimensional resistivity mapping with a computer- controlled array. First Break, vol.8 No. 4, pp.121-129.

Griffiths, D.H. and Baker R.D. 1993. Two-dimensional resistivity imaging and modelling in areas of complex geology. J. Appl. Geophys., 29, pp.211-226.

Hernandez-Stefanoni J.L., Ponce-Hernandez R. 2006. Mapping the spatial variability of plant diversity in a tropical forest: comparison of spatial interpolation methods. Environ Monit Assess. 117(1-3), pp.307-34

Krige, D.G. 1951. A statistical approach to some basic mine valuation problems on the Witwatersrand. J. Chem. Metall. Min. Soc. S.Africa. Vol. 52, No. 6, pp.119-139.

Matheron, G. 1963. Principles of geostatistics. Econ. Geol. 58, pp.1246-1266.

Matheron, G. 1976, Le Choix des Modeles en Geostatistique Advanced Geostatistics in the Mining Industry, NATO Advanced Study Institutes Series Volume 24, pp 11-27.

Minasny, B., McBratney, A.B., and Whelan, B.M., (2002). VESPER version 1.5. Australian Center for Precision Agriculture, McMillan Building A05, The University of Sydney, NSW 2006. (http://www.

Mohmmadi, J. 2006. Pedometrics: Spatial Statistics (Geostatistics). Pelk Publications, 453 p.

Olayinka, A.I. 1999. Advantage of two-dimensional geoelectrical imaging for groundwater prospecting: case study from Ira, southwestern Nigeria. J. Nat. Ass. Hydrogeol., 10, pp. 55-61.

Olayinka, A.I., Yaramanci, U. 1999. Choice of the best model in 2-D geoelectrical imaging: case study from a waste dump site. European J. Environ. Eng. Geophys., 3, pp. 221-244.

Ouadif L., Bahi L., Akhssas A., Baba K. and Menzhi M. 2012. Contribution of Geophysics for the Determination of Aquifers with a Case Study. International Journal of Geosciences, Vol. 3 No. 1, pp. 117-125.

Ouadif, L., Bahi, L. and Baba, K. 2014a. Geophysical Investigations for Groundwater in Outita, Morroco, using ERT Method. International Journal of Engineering and Technology, Vol 5 No 6, pp. 5191-5195.

Ouadif, L., Bahi, L. and Baba, K. 2014b. Identification of formations of soil and subsoil using finite elements modeling. MATEC Web of Conferences pp. 03003-1-03003-4.

Vicente-Serrano S.M., Saz-Sanchez M.A., Cuadrat J.M., 2003. Comparative analysis of interpolation methods in the middle Ebro Valley (Spain): Application to annual precipitation and temperature. Clim. Res. Vol.24, pp.161-180.

Webster, R. 1985. Quantitative spatial analysis of soil in the field. Advances in Soil Science, Vol. 3, pp. 1-70.

Webster, R., Oliver, M.A., 2001. Geostatistics for Environmental Scientists. John Wiley & Sons Ltd, Chichester.

Willmott, C.J. 1982. Some comments on the evaluation of model performance. Bull. Am. Meteorol. Soc.; Vol.63, pp. 1309-1313.

Yans, J., Amaghzaz, M., Bouya, B., Cappetta, H., Iacumin, P., Mouflih, M., Selloum, O., Sen, S., Storme, J.Y, Kocsis, L., and Gheerbrant. E. 2014. First carbon isotope chemostratigraphy of the Ouled Abdoun phosphate Basin, Morocco; implications for dating and evolution of earliest African placental mammals. Gondwana Research 25 pp. 257-269.

Khadija Baba (1), Lahcen Bahi (1), Latifa Ouadif (1)

JGIE Laboratory, Mohammadia engineering School, Mohammed V- Agdal University, Rabat, Morocco)



Manuscript received: 12/12/2013

Accepted for publication: 21/09/2014

Table 1: Summary statistics of random samples

Statistics             AB40         AB80         AB120

number of sample       5151         5151         5151
sum                    808077,22    503542,9     383859,2
mean                   156,877736   97,7563386   74,5212968
Max                    659          400          250
Min                    9,2          7,1          4,1
Variance               10597,0331   2299,88124   1160,37271
standard deviation     102,951886   47,9617328   34,067551
95% confidence limit   0,08995038   0,04190478   0,02975391

Table 2: Summary details of fitted Exponential models

         [C.sub.0]     [C.sub.1]    [A.sub.1]   SDD (1)
        ([m.sup.2])   ([m.sup.2])      (m)

AB40       555,0         13714        33,35     Mod. (2)
AB80       570,0        1991,0        137,3      strong
AB120      444,0         812,0        76,75       Mod.

[C.sub.0] = Nugget effect [C.sub.1] = Sill - Nugget effect
[A.sub.1] = Range
COPYRIGHT 2014 Universidad Nacional de Colombia, Departamento de Geociencias
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2014 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:GEOSTATISTICS
Author:Baba, Khadija; Bahi, Lahcen; Ouadif, Latifa
Publication:Earth Sciences Research Journal
Article Type:Report
Geographic Code:6MORO
Date:Dec 1, 2014
Previous Article:Severity classification of a seismic event based on the magnitude-distance ratio using only one seismological station.
Next Article:Estimating shear wave velocity using acceleration data in Antakya (Turkey).

Terms of use | Privacy policy | Copyright © 2020 Farlex, Inc. | Feedback | For webmasters