1. Introduction
The extraction of gypsum mineral (CaSO4·2H2O) from the Bologna and Romagna Apennines dates centuries (Fabbri et al., 2021). Its use as a construction material can be appreciated in the basements of many historical buildings of the area, among which the iconical “two towers”: Asinelli and Garisenda, the Bologna city symbol (Pesci et al., 2011; Carpinteri et al., 2016; Dallavalle et al., 2022). The specific gypsum type of the area is selenite, with the peculiarity of depositing in layers, and with transparent translucent crystals. The term “selenite” dates back centuries, and it refers to the ancient Greek name of the moon, for its property to reflect moonlight (Del Monte, 2005). Some recent studies proved the resistance of selenite to urban air quality deterioration, with respect to other natural stones (Del Monte et al., 2001).
The use of gypsum at the industrial level in the area dates back to the 1950s: the opening of the “Monte Tondo” quarry between the municipalities of Riolo Terme and Casola Valsenio began in 1958 by ANIC S.p.A., a state company based in Ravenna. The first industrial use of gypsum was dedicated to the production of ammonium sulphate, used in agriculture as fertilizer. Later in the 1970s, a second use as an additive for cement production was found, too. This new application became predominant in the 1980s. Meanwhile, a factory dedicated to the production of premixed plasters, owned by VIC Italiana S.p.A., was realized nearby and notable amounts of the gypsum from the Monte Tondo quarry were sent to the plant. The use of gypsum in the building sector grew rapidly in the 1990s, and the factory, which changed the ownership and passed to the English group BPB, started to produce plasterboards. In 2005, the plant and the quarry were both acquired by the French Group Saint-Gobain. From that moment onwards, all the gypsum extracted from the “Monte Tondo” quarry has been processed internally by Saint-Gobain in the nearby plant, for a variety of building applications: plasterboards, plasters, tiles, sealants, adhesives. All these final products are sold with the commercial brand Gyproc®. The attention on environment and sustainability has posed new challenges to the gypsum valorisation process, and so the life cycle assessment has become a standard in the sector (Pedreno-Rojas et al., 2020; Weimann et al., 2021; Morsali et al., 2024).
With such premises, Saint-Gobain Italia nowadays has the possibility and duty to control all the purity requirements of the raw material, on both the extraction and processing sides. Currently, the excavation method is drill-and-blast, applied over several benches composing the mid-mountain quarry (Casertano et al., 2024). In order to be selected for the plant, the gypsum concentration must be comprised of between 83% and 87% blasted rock. Excavated materials with lower and higher concentration values are stocked in low-grade and high-grade stockpiles, for subsequent mixing. The aim is to provide enough feed for the processing plant at a stable concentration (Merli, 2024).
In this context, it is evident how the interest shifted from merely separating the gypsum from clay levels (which was the common practice) to instead estimating, with the highest degree of accuracy possible, the gypsum concentrations in the different quarry benches and authorized areas for extraction, in order to better control the feed of the plant.
This paper presents the work performed to respond to this need. All the historical borehole data and the recent blasthole data were collected, and the relevant information was extracted and digitalized. The topographic surveys and the authorization limits were integrated to define the intervention area and to create proper boundaries at the stage of 3D modelling. The geostatistical tools of variogram and kriging were then employed to estimate the concentration variabilities within the mining area, together with their accuracy and precision of the estimation, defined by the estimation variance (Matheron, 1971). Such geostatistical tools are commonly applied by the mining industry to natural deposits, but can also support the characterization of stockpiles of excavated raw materials (Kasmaee et al., 2018; Kasmaee et al., 2019).
The methods, even though traditionally applied to large metal mines, have been started to be tested on quarries, for improving the production efficiency (Ozdemir et al., 2018; Busuyi Afeni et al., 2021; Joshi et al., 2022). Moreover, geostatistics can be also applied in other fields, such as to estimate rock properties (Pakdaman & Moosavi, 2023), to map pollutants (Salvi et al., 2023) or to predict natural phenomena (Kasmaee & Tinti, 2018). By applying such geostatistical knowledge, the geological 3D model of the gypsum deposit was then realised using the historical boreholes; it was later combined with the concentration block model, obtained by the information from the blastholes. Geological interpretations and grade variability analyses were subsequentely performed. The final result was the variability assessment of the gypsum concentration in the authorized area, for a better and more efficient planning of the excavations in the upcoming years.
2. Methods
The local geological and geomorphological setting are characterised by over 350 m thick evaporite sequence cut by the Senio stream. The original sedimentary sequence is strongly fractured and has undergone stacking as a consequence of ancient submarine landslides. The geological formation is made of 16 gypsum layers separated by clayey interlayers and deposited during the Messinian salinity crisis, around 6 millions year ago (Lugli et al., 2008). The quarry has a typical geometry of halfway benches and extends from 220 m a.s.l. (base square) to approximately 340 m on the N-E side, intercepting some of the gypsum layers, among which four are currently under excavation by drill-and-blast techniques (see Figure 1). The quarry is composed of fifteen benches of variable height, from 10 to 15 m, and one bench of 20 m height, located at the base square. It is worth noticing that below the base square an abandoned tunnel system is present, heritage of the previous underground mining methods, where bat colonies of various types have settled during the years (Clavel, 2024). It confirms that plants and animal species may thrive in post-mining habitats (Salgueiro et al., 2020). Moreover, besides and below the quarry, some karst cavities of natural importance exist all around the authorized area, for a total length of 6 km, among which the most famous and important is the “Re Tiberio” cave (Sevil-Aguareles et al., 2025). Furthermore, the presence of gypsum in the area favors the natural creation of around sixty sulphur springs, a geological phenomenon occurring in other similar contexts (Vlahovic & Skopljak, 2019). Renowned baths, located a few kilometres down from the “Monte Tondo” extraction site, have been exploiting such waters since 1877. Finally, the quarry is located beside the UNESCO natural heritage site named “karst and evaporite caves of the northern Apennines“. Due to the above-mentioned conditions, the mining activity is strongly conditioned by environmental and nature preservation constraints; therefore, the ongoing exploitation activity, and its expansion, must be managed with extreme caution (see Figure 2).
Figure 1. Eastern view of the quarry, with evidenced the four gypsum layers currently under excavation (Marabini & Poli, 2021)
Figure 2. Satellite image of the quarry (Google Earth - Date: November 2024), encircled by the authorization limits for extraction activities
2.1. Data collection
Various boreholes were drilled throughout the years (1967, 1986, 1994, 1996, 2004) for different purposes, mainly to identify the geological layers and the presence of groundwater. The boreholes of 1967 and 1994 refer to the area of investigation and included the measurement of gypsum concentration; for such reasons, they were used for the present work. Each borehole reported a log with gypsum concentration in the various layers (see Figure 3).
Figure 3. Extract of a borehole sheet from 1994
Lithological data (geological zones) were available within the borehole sheet which were extracted as an input for geological modelling. Based on the input data, four main lithological units were detected: clay, gypsum with small crystals, gypsum with large crystals and sandstone.
Statistical mean was performed to calculate the average gypsum concentration. The boreholes data (altitude, depth, coordinates and calculated average gypsum concentration) are reported in Table 1.
Table 1. General data of the boreholes drilled in the area for sampling, with average gypsum concentration
The localization of the boreholes in the quarry is displayed in Figure 4.
Figure 4. Localization of the boreholes in the quarry
Added to this, starting in 2021, the average concentration of gypsum has been measured on site for each blasthole. During each drilling, a bag of extracted material has been filled directly through the machine. Each bag has then been brought to the quality lab for proper analyses of the properties of interest (see Figure 5).
Figure 5. Sample of results of lab analysis for one blasthole
A sample of the data for the period 2021-2023 is reported in Table 2.
Table 2. A sample of 10 records of recent blasts in the area with average gypsum concentration taken from blastholes
The localization of all blastholes in the years 2021-2023 in the quarry is displayed in Figure 6.
Figure 6. Localization of the blastholes (years 2021-2023) in the quarry
All the information used for the present work were then obtained by combining blastholes and boreholes (see Figure 7).
Figure 7. Localization of both blastholes and boreholes used for the analysis
2.2. Data analysis
As the first step, fundamental for geostatistical analysis, the input data were statistically studied. The statistics parameters of the blastholes are reported in Table 3.
Table 3. Statistics of gypsum concentration and length for all the analysed blastholes
Focusing on gypsum concentration, the percentile analysis is reported in Table 4.
Table 4. Percentile analysis of gypsum concentration in the blastholes
The frequency analysis is reported for gypsum concentration and blastholes length in Figure 8.
Figure 8. Frequency distribution of gypsum concentration in blastholes (left) and their length (right)
As Figure 8 shows, the most representative class of blastholes length is the 16 m class, followed by the 12 m class, and, detached, the 24 m class. The blastholes were therefore divided into three groups, depending on the predominant lengths of the holes: 12 m, 16 m and 24 m. The statistics relating to these groups are reported in the following tables: Table 5, Table 6, Table 7 and Table 8.
Table 5. Statistics of blastholes divided by classes of length
Table 6. Absolute and relative frequency of the gypsum concentration for the class L=12 m
Table 7. Absolute and relative frequency of gypsum concentration for the class L=16 m
Table 8. Absolute and relative frequency of gypsum concentration for the class L=24 m
As shown in Table 5, there are only 27 samples with 24 m length, which can be divided into two 12 m samples and include them in the 12 m class. Therefore, two classes of 12 m and 16 m were considered for the next step analysis.
2.3. Geological modelling
3D geological modelling was constructed using lithological zones, extracted from borehole data in RockWorks® (Rockware, 2024) software. The RockWorks® software was selected being a widespread and robust tool for such types of analyses. The simple method of inverse distance weighting – IWD (power 2) was used for the geological modelling since the lithology was quite homogenous (see Figure 9). The neighbourhood parameters for performing IWD were selected equally to the parameters used for the geostatistical analysis (see Section 2.4).
The 3D geological model included the four main geological zones as:
Figure 9. 3D geological representation (left) and vertical cross-section passing through all boreholes (right) showing the geological model of the gypsum quarry
A unique grid with 5×5×5 m3 dimension was used for both geological modelling and gypsum concentration model. Constructing the geological model is an important step which can help interpreting the concentration block modelling and variation of gypsum concentration based on geological zones’ variation.
2.4. Geostatistical analysis
To use the blastholes, at the center of each blast, two parameters were assigned:
the average gypsum concentration, which is representative of the blastholes grade;
the support (thickness), calculated as the average length of the blastholes, in case they differ from each other; regularization was used to put all blasts as the length of 16 m. The 16-m length was selected based on statistical analysis shown in Figure 8 and from Table 5 to Table 8, and since the highest frequency of data was 16 m. Therefore, few data have been changed (smoothed) by regularization.
In this way, gypsum blast data, each one with cylindrical support, were obtained for the quarry (see Figure 10).
Figure 10. Grades of cylindrical support representing the blasthole groups in the quarry area. The coordinate reference system is local.
The statistics over the cylindrical support are reported in Table 9 and Figure 11.
Table 9. Statistics over the cylindrical support
Figure 11. Histogram of the sixty-three-gypsum grade and thickness data
Variographic analysis was performed over the data. The samples were divided in two groups, those of length 12 m and those of length of 16 m. The corresponding variograms were calculated (see Figure 12).
Figure 12. Experimental variogram, model and number of pairs over the two set of data, those with thickness 12 m (left) and those with thickness 16 m (right).
Due to the data behaviour, number of pairs and field dimension, the same spherical model was chosen to approximate both spatial structures, but using different parameters. The specific parameters used are reported in Table 10.
Table 10. Type and parameters of the variogram models selected
The 16-m variogram corresponds to the spatial structure of the data mode (see Figure 12). Therefore, the corresponding model was chosen for further analysis.
The volumes chosen for the geostatistical analysis were the ones delimited by the Digital Elevation Model (DEM) of the quarry for the years 2022 (DEM2022) and 2023 (DEM2023), and also the final quarry morphology (DEM-final-design). The topography of the quarry is yearly updated, since the company must report every year the excavated volumes to the regional geological and environmental authority.
Ordinary Kriging was performed, using the software ISATIS.NEO® (Geovariances, 2024), to reconstruct the gypsum estimation in two following quarry fields:
The volume comprised between the DEM2022 and the DEM2023, of which the gypsum grade is known. This was used as reference to validate the correctness of the variogram model and the estimation procedure adopted.
The volume comprised between the DEM2023 and DEM-final-design, in order to obtain the final estimation of gypsum concentration in the remaining area to be exploited.
The block model used for the estimation of gypsum concentration has dimension 5×5×5 m3. The details are reported in Table 11, while Figure 13 shows its assonometry.
Table 11. Detail of the grid used in the block model
Figure 13. Block model for the gypsum quarry
3. Results
After proper cutting, the block model for the two quarry limits (DEM-2023 and DEM-final-design) are presented in Figure 14 and Figure 15. The cut between the two is the final block model where performing the estimation.
Figure 14. Block model representing the DEM-2023
Figure 15. Block model representing the DEM-final-design
The estimation of gypsum concentration, together with the associated estimation variance, for the cut block model DEM-final-design – DEM-2023 are reported in Figure 16 and Figure 17.
Figure 16. 3D model of the estimation of gypsum concentration
Figure 17. 3D model of the estimation variance
4. Discussion
Geostatistical 3D modelling is not a common step in small-scale quarries in Italy. Nowadays, based on the trend towards digitalization and the sensitive environmental monitoring, this work was a good opportunity to transform old data and maps into digital versions, which is helpful not only for reporting and social acceptance strategies, but also for re-evaluating the concentration variability in a more precise way. Therefore, the Monte Tondo gypsum quarry is at the forefront in such a context. Within this work, all old data and maps were digitalized and re-evaluated to create a 3D model. There were many challenges to create such 3D models, since most of quarry process did not consider important points for estimation, such as support of sampling, spread borehole grid, etc.
The estimation map of the quarry (see Figure 16) shows a major gypsum concentration in the South West benches, which is coherent with the four gypsum layers represented in Figure 1. On the other hand, less concentration values of gypsum are predicted on the North East benches, where large volume of clay at levels 290-310 m have been effectively found. Such a presence of clay is coherent with the geological model of Figure 9, where clay was defined at the top layers of the model. The estimation variance map presented a low value in the middle of the block model, while higher values can be found at the model borders. It highlights a decrease of precision of the estimation from the centre to the borders.
By comparing the model of the estimation variance (see Figure 17) with the 3D view of samples (see Figure 10), it is clear the precision is greater near the samples and where information is available, i.e. mainly where the blasts were carried out between 2021 and 2023. The high density of the data has made it possible to reconstruct a reliable representation of the gypsum content in the exploitation areas of the quarries.
Concerning the volume still to be excavated (corresponding to 372,100 m3, comprised between DEM-final-design and DEM-2023), based on the carried out estimation, it is expected to get an average gypsum concentration around 87.2%, with an average estimation standard deviation around 0.7%. Such results indicate high theoretical precision of estimation. The gypsum concentration, in average, is at the upper border of the correct range for feeding the processing plant. The highest values are expected in the southern area of the quarry and at greater depths, while the lowest values can be found in the northern area and at higher elevations. This means that almost all gypsum benches can be proficiently exploited, in order to realize a correct and stable concentration mix to the processing plant. The result is highly beneficial, since it provides evidence to the regional authority of the potential and the goodness of the quarrying activity, with minimum stock of mining residues in the environment. The details of gypsum concentration forecast for the different benches are reported in Table 12.
Table 12. Average values of the estimation, for each quarry bench
The comparison between the geological model and the concentration block model was helpful to interpret the variation of gypsum concentration in the different zones and benches of the quarry. The low concentration zones of gypsum (evidenced by the estimated block model) are linked to the clay zones (represented in the geological model). This comparison confirms the variability of gypsum concentration, which is consistent with the geological changes.
Thanks to the estimated block model, the company has the possibility to further develop the production schedule of the quarry, and to subsquentely optimize the excavation and crushing process, in order to get the correct and stable gypsum concentrations, with minimum loss of raw materials. Such activities are an international mining standard, and it is worthy and promising they start to be applied in quarry operations, too.
5. Conclusions
Digitalization is a fundamental stage for all mines and quarries in Europe and one of the strategical steps that can help even small-scale quarries to have a better understanding of their geology and concentration models. This is possible by inserting all historical data (boreholes) and active data (blastholes) in consistent softwares, in order to create robust geological and concentration models. The quarry operations in the “Monte Tondo” case study benefited from digitalization and geostatistical modelling to construct both geological and gypsum concentration models. The realization of a detailed 3D model representing the estimation of the gypsum content in the target exploitation area allowed for the improvement of the quarry management. This model is obtained by blastholes, while boreholes contributed to creating the geological setting, which helped in interpreting the lithological changes. Specifically, clearly distinguishing the low concentration areas from the high ones allowed for the efficient management of the raw materials mix, as to keep the inlet of the processing plant always inside the admissible concentration range.
Moreover, the 3D visualization allowed for better monitoring of the progress of the excavation, how the morphology of the quarry evolved over time and then to tackle, sufficiently in advance, any type of geological situation different from what was originally expected.
Author’s contribution
Stefano Merli (MSc): data curation, investigation and visualization. Sara Kasmaee (PhD): conceptualization, data curation, formal analysis, methodology, validation, visualization, writing – original draft and writing – review & editing. Cristina Barbalucca (MSc): data curation, formal analysis, resources and supervision. Vanessa Cellini (MSc): data curation and project administration. Francesco Tinti (PhD) conceptualization, methodology, project administration, writing – original draft and writing – review & editing.
All authors have read and agreed to the published version of the manuscript.