Phosphorites appear in three principal zones of Tunisia: the Northern Basins, the Eastern Basins (NortheSouth Axis area), and the Gafsa-Metlaoui Basin. The Phosphorites are mainly present in the Paleocene-Eocene chouabine formation and its lateral equivalents within the Metlaoui group; they form part of the Middle Eastern to North African Late Cretaceous-Paleogene phosphogenic province     .
The economic deposition of phosphorite is mainly in the Gafsa-Metlaoui basin in central Tunisia. Here, phosphorites are operated by the Compagnie des Phosphate de Gafsa (CPG) and are considerably utilized as raw materials for fertilizer production by the Groupe Chimique Tunisien (GCT). In other Tunisian basins, phosphorite layers are usually thinner and show inferior quality. Anterior studies have discussed stratigraphy, mineralogical composition, sedimentary facies, deposition environment and the diagenetic history of Tunisian phosphorites  -  , but reserve estimation of phosphate remain limited. Published data are predominantly by geochemical analysis of the Gafsa-Metlaoui deposits.
An accurate estimation of the ore reserves of a specific mineral is required in order to better manage and plan its extraction at mines. A thorough understanding of the distribution of the mineral content of the ore is also necessary. Nowadays, with the help of different types of computer software and through the use of geostatistical techniques, estimations can be as accurate as required. However, the implementation of these modern tools has to be done logically and with sufficient data.
Phosphate ore is known to be one of the chief minerals produced in Tunisia. The Tozeur-Nefta deposit is found in the southwestern part of the Gafsa mining basin which is approximately 12 km west of Tozeur. Phosphate rocks belong to the extensive Paleocene-Eocene phosphate series and are widely used in the chemical industry, as a source of phosphor in fertilizers and in manufacturing chemicals such as phosphoric acid . The different grades of phosphate rocks can be distinguished according to the following: Poor-quality grade, P2O5% which ranges 15 to 22; medium grade, P2O5% ranging from 22 to 27; and rich grade, P2O5% greater than 27 . Thus, identifying the geographical distribution of phosphate is important for engineers in this mining industry and would assist in choosing the most appropriate mining method for the extraction and for production control.
2. Geographical Setting
The main feature of Tunisia is that it has two varied geological regions: the north is characterized by the folded and faulted Atlas Mountains while the south, the stable Saharan platform     . The area explored in this study, known as Draa El Jerid (the Nafta-Tozeur phosphate deposit) is located in the Gafsa-Metlaoui Basin (Figure 1).
This basin is situated in the southern Atlas of central Tunisia, with an area of around 4500 km2 . From the structural point of view, it is a transitional region separating an actively folded and faulted zone to the north, the central-northern Atlasic Domain, and the stable Saharan Platform to the south   . The Basin is limited in the north by the Metlaoui mountain chain, containing Jebels Ben Younes, Bouramli and Orbata, and in the south by the North Chotts Range.
Figure 1. (a) Location map of the study area; (b) Geological map of the Tozeur-Nefta deposit .
In this area, a noticeable subsidence has created two hypersaline depressions called El Gharsa Chott, El Jerid Chott parted by the northern uplands of the Chott and the Draa Jerid anticline also known as the Tozeur uplift or Tozeur Ridge    (Figure 2).
In the Gafsa-Metlaoui Basin, the anticline structures of Draa Jerid, located midway between El Gharsa Chott and El Jerid Chott (Figure 1), make up a portion of the western part of the Chotts fold belt   , corresponding to structures in the most southerly part of the Atlassic domain     .
Through the study of geoseismic cross section X Guellala et al.   concluded that in the western edge of the Nafta area, there is a distinct discrepancy in the depth of the lithostratigraphic formations between the geological structures. Within the anticline of Draa Jerid, the Sidi Aich and Boudinar formations are respectively located at 2100 m and at 2430 m. These formations are respectively reached at a depth of 2850 m and at a depth of 3150 m within the Jerid basin.
In addition, this cross section reveals a blocked communication between Sidi Aich and Boudinar aquifers. This is explained by the presence of the East-West reverse fault in the North and in the South of the Draa Jerid anticline (Figure 3(a)).
The most noticeable difference between the piezometric values of the Continental Intercalaire aquifer detected at Hezoua and those measured in Nafta  (Figure 3(b)), reflect a blocked groundwater circulation between the two regions, expression of a normal fault which subsides the northeastern compartment and puts the Sidi Aich formation in the anticline of Hezoua at the front of Orbata and Zebbag formations in the western boundary of the Draa Jerid (Figure 3(b)).
At the east of the Draa Jerid area, based on MR borehole, TZ borehole and the geoseismic cross section V, the Sidi Aich aquifer formation, reached at a depth of 1735 m within the TZ borehole is laterally in front of the Zebbag formation of El Jerid Chott and in front of the Orbata and Zebbag formation of El Jerid basin
Figure 2. Cross-section of Tozeur uplift in the Jerid basin .
Figure 3. (a) Geoseismic cross section X western extremity of Draa Jerid  ; (b) Normal fault: western extremity of Draa Jerid; (c) Geoseismic cross section V eastern extremity of Draa Jerid .
. In El Jerid bassin, the MR borehole reaches the Sidi Aich formation at a depth of more than 2500 m. This is explained by the presence of the east-west reverse fault in the north and in the south of the Draa Jerid anticline (Figure 3(c)).
The oriental extremity of Draa Jerid is marked by the NW-SE fault of Negrine-Tozeur, the western branch of the South atlasic accident in Tunisia  . This fault extends from the Negrine to Kebili area via Tozeur  west of the Jebel Sidi Bouhlel megastructure .
Therefore, the Nefta-Tozeur area appears as an isolated structure, limited by two reverse E-W faults on the northern and southern sides and two normal NW-SE faults on the eastern and western sides (Figure 4).
The Gafsa-Metlaoui Basin is mainly range in age from Cretaceous to Quaternary. Deposition in the Gafsa-Metlaoui Basin developed in a partially confined environment, which oscillated from littoral to lagoonal atmosphere, conducing to regular or periodic deposition  . As a result, there is notable diversity of facies inside the basin, containing porcelanites, phosphorites, shales, cherts, limestones, marls, gypsum and dolostones . Marly intercalations inside the phosphorite series usually comprise from 1% to 2% until 7% of TOC (Total Organic Carbon). The organic matter is immature and is mainly of bacterial origin and marine planktonic   .
The principal phosphorite sequence composes the Chouabine member   of the Metlaoui Group of the Paleocene-Eocene . The Group reposes
Figure 4. Simplified structural map of Draa Jerid structure.
over the Maastrichtian-Danian El Haria Formation, and is covered by Jebs Formation of the middle to upper Eocene   .
The Metlaoui Group is divided into 3 members:
• A lower carbonate and evaporitic member, dated from the Palaeocene (Selja member);
• The main middle phosphatic Chouabine member, attributed to the Upper Palaeocene, is represented by ten phosphatic deposits interbedded with carbonates, clays and a cherty level;
• The upper carbonate member, named Kef Eddour, is comprised of bioclastic limestones, covered by a phosphatic recurrence named “the Upper Phosphate”, and is overlain by limestones and dolomites.
The Chouabine member regularly reaches 25 to 100 m thick    . Ten essential phosphorite layers (from top to bottom: layers 0-IX) constitute this member, isolated by layer of marly limestone, marl, chert and diatom-bearing porcelanite . In the eastern part of the basin the quantity of phosphorite diminishes and an important amount of carbonates and marl develops  , however to the west biosiliceous deposits (diatom-bearing porcelanite) become more significant    . Diatom faunas mark a warm climate, coastal shallow-marine context within the central basin, with more salty eutrophic circumstances to the east . Water depths augmented for east to west  , reaching a maximum water depth of 100 m . The phosphorite-organic-rich marl-diatom-bearing porcelanite facies indicates the classic coastal upwelling trinity . The oyster-rich limestones with phosphorite interbedding constitute the upper part of the Metlaoui Group, recognized by the miner as the “phosphate du toit”; these layers represent the Kef Eddour Formation .
The Phosphorites of the Gafsa-Metlaoui Basin were examined by Garnit  at five localities which are the Naguess deposit at the north of Jebel Alima, the central Kef Eddour deposit which is situated around 10 km NNW of Metlaoui, the Table Metlaoui deposit is situated on the southern side of the line of hills ranging between Jebel Stah and Jebel Alima (Figure 1). All these sites show the average total thickness of the phosphorite beds between 7.2 m and 12.4 m, with about average total thickness of the intercalated layers between 11.3 and 12.8 m.
Precise dating of the Chouabine and analogous phosphorites has confirmed challenging. Particular chemostratigraphic and biostratigraphic analyses have been studied just on Gafsa-Metlaoui Basin deposits. In the main, the begining of phosphorite sedimentation has been judged as a marker of base of Ypresian .
Ben Abdessalem  located the Paleocene-Eocene limit at the base of layer II in the Gafsa Metlaoui series, relying on the appearance of an organic walled dinoflagellate cyst aggregation predominating by Apectodinium spp. in layer 0-I.
Bolle et al.  fixed the base of the Eocene at the summit of the Chouabine Formation in the GafsaMetlaoui Basin, on the basis of an experimental assignment of layers 0-II to calcareous nannofossil zone NP9 with records of Discoaster multiradiatus Bramlette and Riedel, and a sequence stratigraphic correlation to the Elles section in northern Tunisia.
Ounis et al.  and Kocsis et al.   locating the base of the Eocene at the summit of Chouabine Formation in layer III in the Gafsa-Metlaoui Basin.
Some time ago, El Ayachi et al.  inserted the Paleocene-Eocene limit in the inferior layers of the Chouabine Formation in the Oued Thelja section, relying on two samples that yielded planktonic foraminifera.
The Tozeur-Nafta phosphate deposit, found in Draa Jerid, is derived from the phosphatic Gafsa-Metlaoui basin and this deposit is alike other deposits in the Gafsa phosphate basin.
Within the Tozeur-Nafta deposit sixty two (62) boreholes have been sampled (Figure 5), with depths ranging between 30 and 1100 m, to investigate more thoroughly the phosphate deposit of the Metlaoui Formation and to attempt to estimate the reserves present in this location.
Figure 5. Borehole distributions in Tozeur-Nefta deposit.
Figure 6. (a) E-W deep wells correlation in Tozeur-Nefta deposit; (b) 3D stratigraphic model of Tozeur-Nefta deposit; (c) Well location map; (d) E-W cross section of Tozeur-Nefta deposit; (e) N-S deep wells correlation in Tozeur-Nefta deposit.
within the Metlaoui Formation, only the middle phosphatic Chouabine and Kef Eddour members appear covered by the Beglia formation attributed to the Miocene.
This correlation also indicates that, in the Tozeur-Nafta deposit from base to top, CIII, CIV and CV do not appear. In this correlation, the phosphatic Chouabine member shows a marly layer at the bottom with 8.50 m of thickness surmounted by grey to a brownish phospharenite layer (CIX) which make non-productive layers because of the high content of MgO (1.79%). The centimetric to decimetric thick layers (CIII, CIV and CV), are considered sterile owing to their slight thickness and six main metric economic phosphorite beds are exploited (C0, CI, CII, CVI, CVII and CVIII) separated by marly, carbonated and chert levels. The total thickness of the phosphorite beds averages 9.8 m and the intercalated beds, 14.95 m (Table 1 and Table 2).
Regarding the thickness of the afore-mentioned, the three main phosphate layers are CI, CII and CVI, which contribute to 74.5% of the total thickness (7.30 m). In general, the phosphate layer thickness in the northeastern limb of the deposit S.49 (8.71) is more evident than in the southeastern limb S.70 (6 m) (Figure 6(e)).
The maximum thickness of the phosphate layers is situated on those in the eastern part of the deposit (Figure 6(b) and Figure 6(d)), in particular around S.9 (12.22 m) and S.30 (11.45 m) boreholes with the exception of the S.69 (8.16 m), S.68 (8.95 m) and S.28 (8.82 m) boreholes.
Concerning the Kef Eddour carbonates formation (Upper member of the Metlaoui group), the upper unit is marked by two thick carbonate bars comprised of colossal amounts of bioclastic dolomitic limestone with phosphorite intercalations, recognised in this mining industry as the “phosphate du toit” . The Kef Eddour Formation show an average thickness of 40 m. The thickness is more noticeable towards the west (Figure 6(b) and Figure 6(c)).
The Beglia formation (Mio-Pliocene) mostly involves sandy fluvial deposits. The thickness of the afore-mentioned ranges from 0.00 m (P.2) and 3 m (S.25) to 105.00 m (S.51), showing an average of 21.00 m and is generally thicker to the west (Figure 6(b) and Figure 6(d)).
To generate a continuous surface, it necessitates vectorial data (points, lines, surfaces) containing value information, which is interpolated in order to ensure a continuous surface. The precision of the model depends on the used methods of interpolation. Therefore, it is essential to examine, the performance of each technique.
The first stage of this was to compare the general techniques of interpolation of: Kriging, IDW and Spline, keeping in mind the knowledge that different interpolation methods show varying strengths and weaknesses depending on the
Table 1. Average thickness of exploitable phosphate layers.
Table 2. Average thickness of the interlayers.
dataset used. There can be no generalization made on whether one interpolation method would perform better than another without taking into cognition, the type and nature of the dataset and Implicated phenomenon.
The aim of this comparison was to experimentally assess the performance of Kriging, IDW and Spline methods of interpolation in approximating unidentified phosphate thickness values and modeling spatial distribution of phosphate deposit. Through this experimental evaluation, a comparative analysis was made based on the prediction mean error, prediction root mean square error (Figure 7) and validation outputs of these methods of interpolation. The investigative results for each method on both biased and normalized data show that the most accurate interpolation in this experiment was provided by the Kriging technique.
Kriging, described as “interpolation with geostatistics”  , is a technique used to analyze continuous data. Interpolation using this technique is widely approved by statisticians and scientists since the process is based on theory. There are various kriging techniques available through different software such as GSLIB, Geostat, ArcGIS, ArcInfo, Surfer, etc. The phosphate reserves at the 62 boreholes (Figure 5) have been spatially studied through the thorough characterization of the autocorrelation and semivariogram mechanisms via geostatistical implements.
In this research autocorrelations and correlograms were calculated and constructed using ArcGIS which explored the spatial distribution of phosphate on a point scale system by means of the ArcGIS-Geostatistical analyst.
Figure 7. Prediction errors of the phosphate thickness estimation to evaluate the performance of interpolation methods (a) Kriging, (b) IDW, (c) Spline interpolation methods.
In using Kriging, it was possible to complete the interpolation or appraisal of values for points in an area that had not yet been sampled (this is known as a nearest-neighbor technique). In this case, the values at locations close to the interpolation point are used to estimate the point value of the interpolation. In a data set of the 15 closest neighbours, Kriging uses the core spatial relationship that exists among them. Kriging was applied in this based on a regionalized variable theory and it appears more better than any other interpolation method since it provides the most favoured interpolation estimate for a specific coordinate location, in addition to a variance estimate for the interpolation value.
In utilizing this method, 2D maps of the phosphate content data were created to display the horizontal spatial distribution of phosphate in each depth interval.
Finally, these kriging maps were availed of to estimate the phosphate reserve using methods considered to be advanced statistically and geostatistically. Geostatistical Analyst (version 10.2) provides well-established geoprocessing tools and serves as an extension for advanced surface modeling by means of deterministic (non-geostatistical method) and geostatistical methods.
In the current study, the reserves were estimated with consideration to the thickness of the phosphate layers in the boreholes. The reserve was valued using the geostatistical method as Ordinary Kriging Method (OKM).
At present, there are several computer programs that allow, through geostatistical analysis, the estimation of ore reserves to a high precision, in less time and at a lower cost than others. If the ArcGIS software is used, the interpolation is usually more precise as it relies on the geographical dispersal of the information (actual spatial coordination). In this study, before to using the interpolation techniques, the data was examined using the exploratory spatial data analysis tools alongside the ArcGIS software.
In this case study, sixty two borehole samples were gathered from different locations in the Tozeur-Nafta region so as to analyze and map the spatial distribution of the phosphate thickness. The phosphate thickness kriged maps as seen in (Figures 8-10) indicate that the eastern part of the area under study and the small region towards the west have a greater thickness than the other areas which have a relatively lower thickness. The lowest thickness is shown in yellow.
The phosphate reserve estimation was based on the Equation (1) (Formula source is the Compagnie des Phosphate de Gafsa (CPG)):
where BTS is a raw phosphate, sorted and dried (in tons), VPP is a volume of phosphate set up (m3), VPE is a volume of loss on exploitation (m3), VS is a volume of fouling (m3), VH is a volume of humidity (m3).
VPP calculation was based on the Equation (2):
VPP = Kriging Mean value (m) × Area (m2) (2)
where the Kriging mean value is the mean thickness of phosphate layers between
Figure 8. (a) Predicted thickness map of phosphatic layer C0 in Tozeur-Nefta deposit; (b) Predicted thickness map of phosphatic layer CI in Tozeur-Nefta deposit.
two consecutive isopack curves (m) estimated from the geospatial analysis (Figures 8-10). The area was evaluated from the ArcGIS software. VPE calculation was based on the Equation (3):
VPE = Sr × Tr (3)
Tr is a recovery rate is considered at a loss of 5 cm from the wall and 5 cm from the roof of each layer.
VS calculation was based on the Equation (4):
VS = Area × Ts (4)
Ts is a fouling rate is dependent on the approximated thickness of limestone blocks or the existence of Marl within each layer.
Figure 9. (a) Predicted thickness map of phosphatic layer CII in Tozeur-Nefta deposit; (b) Predicted thickness map of phosphatic layer CVI in Tozeur-Nefta deposit.
VH calculation was based on the Equation (5):
VH = VPP × Th (5)
Th is the humidity equals 12%.
d is the density of the level of phosphate (Table 3).
Tables 4-10 summarize the calculations made for each phosphate layer in order to obtain the phosphate volume from each layer as well as the total phosphate volume in the whole deposit.
The estimation of the reserve of a deposit, based on borehole surveys, leads to the calculation of the ore it contains. These values correspond to the maximum
Figure 10. (a) Predicted thickness map of phosphatic layer CVII in Tozeur-Nefta deposit; (b) Predicted thickness map of phosphatic layer CVIII in Tozeur-Nefta deposit.
Table 3. Phosphate density by layer.
Table 4. Phosphate ore reserve estimations using the Kriging method for phosphatic layer 0 in the Tozeur-Nefta phosphate ore deposit.
Table 5. Phosphate ore reserve estimations using the Kriging method for phosphatic layer I in the Tozeur-Nefta phosphate ore deposit.
Table 6. Phosphate ore reserve estimations using the Kriging method for phosphatic layer II in the Tozeur-Nefta phosphate ore deposit.
Table 7. Phosphate ore reserve estimations using the Kriging method for phosphatic layer VI in the Tozeur-Nefta phosphate ore deposit.
Table 8. Phosphate ore reserve estimations using the Kriging method for phosphatic layer VII in the Tozeur-Nefta phosphate ore deposit.
Table 9. Phosphate ore reserve estimations using the Kriging method for phosphatic layer VIII in the Tozeur-Nefta phosphate ore deposit.
Table 10. Phosphate ore reserve estimations using the Kriging method for all phosphatic layers in the Tozeur-Nefta phosphate ore deposit
potential of the deposit: which are the resources. It should be noted the result of such a calculation is called inferred reserves. After this initial calculation and based on the industrial applications and characteristics of the ore, there can then be a calculation of the recoverable reserves undertaken.
The creation of the different iso-thickness maps has allowed us to conduct an analysis in terms of thickness but also in terms of the financial viability of the different exploitable phosphatic layers. The present study has allowed us to deduce that layers I and II of the deposit are economically very profitable, given that, for example, layer I has a significant thickness can reach up to 7.5 m. The variation of the thickness could be due to the paleogeography of the basin, as well as to the variations of the eustatic sea level during the Paleocene-Eocene. The iso-thickness maps have allowed us to identify the significance of the amount of phosphate found in each analysed layer. Finally, the study of the reserve estimate shows that layers I and II have a reserve of 235 million tons of phosphate, which means that these layers are economically exploitable.
The Paleocene-Eocene phosphate series constitutes the continuity and the equivalent series in the Gafsa mining basin which is situated underground. The thicknesses of the upper layers average about 78 m and consist essentially of sands from the Beglia Formation and carbonates from the Metlaoui Formation. This phosphate series is characterized alternately by phosphate levels and marl levels .
It is known that with heavy investments and fluctuating markets, a mining group needs precision in order to make decisions. Geostatistics has forged a practical tool that enables one to go as rationally as possible from some numerical information to ore “benefited” by crossing the smoothing effect of the geometry of mining operations.
ArcGIS software provides mining companies with the geographic advantage to target mineral potential. Geoscientists of mineral exploration use various types of datasets to search for deposits; such information needs to be viewed and analysed quickly and easily. ArcGIS gathers this data using easy-to-use software applications and tools, making the spatial context of the information available to mining planners and enabling a more in-depth understanding of the geography of potential sites. Using ArcGIS geostatic tools, it is possible to locate profitable deposit regions and calculate the potential of mineral sites through borehole data.
Methodology integrating kriging and GIS has been introduced and implemented in the Tozeur-Nafta area to demonstrate its usefulness in mapping the spatial distribution of phosphate thickness in a specific mining area. In the present study, ordinary kriging was used to identify the different areas consisting of economical mining grades of phosphate and poorer quality ones. Furthermore, GIS was used to be illustrative of the results and to provide more details to the kriging results. From the prediction maps of phosphate thickness, the phosphate quantity in the Tozeur Nefta deposit can be calculated at an average of 325 MT.
From the prediction maps of phosphate thickness, it can also be concluded that, the thickness generally increases in the E-W direction, signifying that the mining engineer would probably start the extraction operations in this economical part.