Automatic Characterization of Block-In-Matrix Rock Outcrops through Segmentation Algorithms and Its Application to an Archaeo-Mining Case Study

Abstract

The mechanical behavior of block-in-matrix materials is heavily dependent on their block content. This parameter is in most cases obtained through visual analyses of the ground through digital imagery, which provides the areal block proportion (ABP) of the area analyzed. Nowadays, computer vision models have the capability to extract knowledge from the information stored in these images. In this research, we analyze and compare classical feature-detection algorithms with state-of-the-art models for the automatic calculation of the ABP parameter in images from surface and underground outcrops. The outcomes of this analysis result in the development of a framework for ABP calculation based on the Segment Anything Model (SAM), which is capable of performing this task at a human level when compared with the results of 32 experts in the field. Consequently, this model can help reduce human bias in the estimation of mechanical properties of block-in-matrix materials as well as contain underground technical problems due to mischaracterization of rock block quantities and dimensions. The methodology used to obtain the ABP at different outcrops is combined with estimates of the rock matrix properties and other characterization techniques to mechanically characterize the block-in-matrix materials. The combination of all these techniques has been applied to analyze, understand and try, for the first time, to model Roman gold-mining strategies in an archaeological site in NW Spain. This mining method is explained through a 2D finite-element method numerical model.

1. Introduction

A great number of geological formations are structurally complex and heterogeneous, often requiring simplifications for their mechanical characterization. A paradigmatic example of this complexity is exhibited by a specific group of materials composed of blocks, which sometimes exhibit relevant variations in shape and size, embedded in a finer-grained, mostly homogeneous matrix.

In nature, these “block-in-matrix” geomaterials can be found in a wide range of lithological units and scenarios, in the form of colluvium and alluvium deposits (Figure 1a) and breccias and conglomerates (Figure 1b), but also as decomposed granitic rock masses (DGs) (Figure 1c) and volcanic agglomerates, pyroclasts, olistostromes and melánges (Figure 1d,e). To designate such materials, Medley [1] introduced the non-geological word bimrock, an acronym standing for Raymond’s [2] block-in-matrix rock, to indicate “a mixture of rocks, composed of geotechnically significant blocks within a bonded matrix of finer texture”. Later, following Medley’s line, the acronym bimsoil (block-in-matrix soil) was used to designate geological units with rock blocks embedded in a soil-like matrix [3,4,5].

Several authors have carried out laboratory and in situ tests, as well as numerical analyses, and have demonstrated that bimunits (from block-in-matrix units [6]) are strongly affected by the quantity, position, shape, orientation and dimension of blocks [7,8,9,10,11,12]. Specifically, when the block content falls in the range of about 20% to 75%, the strength increases, the deformability decreases and failure surfaces become more tortuous.

Due to the aforementioned reasons, a rigorous estimate of block content through the volumetric block proportion (VBP) parameter is essential for the mechanical characterization and modeling of a block-in-matrix material. Moreover, when excavating these complex geomaterials, the presence of unexpected rock blocks can induce, among several other effects, damage to cutters, face instabilities and obstructions [13]. At a rock mass scale, the extraction of representative samples is virtually impossible from a practical perspective. That is why some methodologies have been developed relying on fundamental stereological principles [4,14,15], which state that, when sufficient data are available, VBP estimates can be obtained through 2D (surface), 1D (linear) or even 0D (point) measurement approaches [16,17,18,19].

These dimensions are associated with standard characterization techniques, like imagery, UAV-based photogrammetry [20,21] and geological mapping, drilling and scanlines, and grid-based sampling, and consist in segmenting the dimensional space by counting and measuring elements within it. This task can be particularly difficult for granular materials with large grain-size variability, as stated by Wolcott and Church [22], in sediments with a fluvial origin, such as those analyzed in this paper.

Computer vision techniques, especially focused on segmentation and image analysis approaches, have long since emerged as complementary characterization tools in multiple fields of geoengineering [23,24,25,26] and have proved to be particularly effective in the characterization and sizing of granular geological samples.

Considering the foregoing, the main objective of this paper was to check the capabilities of conventional edge- and corner-detection methods in determining areal block proportions (ABPs) from photos of rock walls excavated in block-in-matrix materials, both in surface and underground openings. Alternatively, a state-of-the-art segmentation algorithm (the Segment Anything Model, SAM) [27] was explored and found to effectively estimate ABPs from several photos showing different features related to the texture, the lighting conditions and the contrast. The algorithm was implemented in a Python v. 3.10 programming language script. The tool was enhanced with the ability to generate histograms representative of the relative block size distribution for the analyzed photo. The results given by the method after the analysis of 10 selected images were contrasted with several estimates visually carried out by a group of experts, a practice followed in other geomechanics-related studies [28,29].

Figure 1. Examples of block-in-matrix materials: (a) detailed view of a typical geomaterial found in a slope in Las Médulas; (b) conglomerate boulder in Covarrubias, Spain; (c) decomposed granite (DG) outcrop in Galicia, Spain, where big boulders are mixed with soil material; (d) tectonic mélange material in a core sample; and (e) a tunneling operation through a block-in-matrix structure with dug-out blocks. Photo credits: (a) taken by the authors; (b) courtesy of Alejano et al. [30]; (c) courtesy of Alejano and Carranza-Torres [31]; and (d,e) adapted from Adam et al. [32].

Figure 1. Examples of block-in-matrix materials: (a) detailed view of a typical geomaterial found in a slope in Las Médulas; (b) conglomerate boulder in Covarrubias, Spain; (c) decomposed granite (DG) outcrop in Galicia, Spain, where big boulders are mixed with soil material; (d) tectonic mélange material in a core sample; and (e) a tunneling operation through a block-in-matrix structure with dug-out blocks. Photo credits: (a) taken by the authors; (b) courtesy of Alejano et al. [30]; (c) courtesy of Alejano and Carranza-Torres [31]; and (d,e) adapted from Adam et al. [32].

Given the potential performance and ease of the algorithm in determining ABP values from a large set of inputs, this work attempts also to provide assistance in future improvements of 3D (VBP) estimates in block-in-matrix materials from images, a topic still deserving further research.

2. Geological, Geomechanical and Computational Framework

In this section, we briefly contextualize the archaeological site under study and present the image database utilized in this research work and the camera used to capture the images. Additionally, the two approaches used to obtain the ABP of each picture, corner–edge detection and image segmentation, are also introduced, together with the SAM algorithm and its implementation.

2.1. The Archaeological Site under Study

Las Médulas (León, Spain) is thought to be one of the most important gold exploitations of the Roman Empire (Figure 2), covering an area of 600 ha and producing 4–5 t of gold from mined material [33], nearly 1 Mm³, which formed tailings that also cover an additional 600 ha of space. The mining technique used involved a complex hydraulic network that required more than 700 km of canals, including local tunnels.

2.2. Geological and Mining Features

Geologically, the area is formed of a sedimentary series associated with the formation of alluvial fans resulting from fluvial deposits from the surrounding mountain ranges [34] and deposited in an environment with an initially warm–humid seasonal climate, with increasing aridity over time. Among these deposits, three main conglomerate formations—resting discordantly on a Paleozoic series [35], corresponding to middle facies of alluvial fans formed in a depositional basin, probably from the Miocene epoch—can be found.

These conglomerate formations (Figure 3 and Figure 4) correspond, from bottom to top, first to the Orellán Formation, made up of alluvial and colluvial deposits containing cobbles up to 10 cm in diameter, embedded in a relatively soft matrix. The thickness of this formation may vary from 20 m to practically no outcropping. The Santalla Formation appears on top of Orellán and is also composed of conglomerates embedded in a quartz, muscovite and clay matrix, reaching a maximum thickness of approximately 25 m.

The maximum size of the blocks, predominantly made of quartzite, is about 1 m in diameter (Figure 4). This formation contains the gold sought by the Romans in quantities of up to 1 g/m³. The Las Médulas Formation lies on top of the previous ones, and its thickness is about 100 m, consisting of alternating beds of boulders and gravel, fine sand, and silt. This material appears to be more resistant, embedded with large blocks, showing a typical deposition in alluvial fans in a high-energy environment. Nevertheless, some finer-grained layers may have appeared as a result of possibly calmer geological events. The Au grade in the Las Médulas Formation is much lower than in the underlying: 10–20 mg/m³.

The development presented in this paper arose as part of a combined engineering–archaeological study to better understand the way in which the Romans mined the gold. For these tasks, the thick cover of the Las Médulas Formation (considered by the Romans as sterile material) had to be removed in order to reach the higher Au grade of the Santalla Formation. In the 1st century AD, Pliny “the Elder” named the technique to remove this material with the illustrative name of Ruina Montium: the “collapse of the mountains” [37]. It indeed was described as the caving of large amounts of material to remove the upper levels and reach the lower ones with a higher gold grade. To investigate the methods employed by the Romans, a geomechanical characterization of these block-in-matrix materials was conducted, relying on a reasonably accurate estimation of the VBP of the conglomerates, something that was eventually performed for different outcrops of the site from ABP estimates, according to the approach herein presented.

2.3. Methods for Bimunit Characterization

Both in the Las Médulas and Santalla Formations, large boulders can appear. Additionally, the matrix is made up of fine materials (fine sand, sometimes silty) that, apparently, are not very strong. All these aspects, in addition to the heritage protection status of the area, prevent taking representative large samples of these materials so that they can be tested in the laboratory or in situ. Therefore, characterizing the behavior at a real scale (such as how rocks respond to excavations and slopes) is a challenging task. Nevertheless, there are not many relevant studies that allow rigorous characterization of these materials. Three types of approaches have been proposed to characterize the properties of these bimunits, which are hereinafter described.

The first type would imply experimental studies on samples [38], in case they could be carved or created artificially, something that does not occur in the case under scrutiny. The second type, used in some large engineering works, would imply carrying out large-scale in situ direct shear tests [9,39], so it would be difficult to implement this in Las Médulas, both in economic terms and because of the damage it would cause. A third, more general type [4,5,8,40] would be the characterization of the matrix and the quantity and nature of rock blocks and, through empirical or numerical approaches, the estimation of the properties indirectly based on the VBP. Particularly suitable for a first characterization approach is the method proposed by Sonmez et al. [40], updated by Kalender et al. [4], suggesting an estimate of bimrock strength based on the rock matrix characterization, VBP, the angle of repose of the boulders, and estimates of adhesion between the boulders and the matrix.

2.4. Camera and Photo Database

Several pictures were taken at different representative locations of the site, using an iPhone 11 smartphone model. These pictures have a size of 1536 × 2048 pixels and a resolution of 72 dpi in JPEG format. Although the quality and resolution of a smartphone camera are lower than those of a professional one, this type of integrated camera has been successfully used by geopractitioners for multiple applications [41,42,43]. Moreover, the obtained pictures are of sufficient resolution, allowing a detailed observation of each scenario. For the present study, 10 photos of different rock slopes taken at several areas of the site were selected.

2.5. Corner- and Edge-Detection and Watershed Algorithms

The initial approach considered for determining the area represented by the blocks in each photograph involved the use of corner- and edge-detection algorithms, applied to the selected images as described in Section 2.3. Complementarily, the capabilities of a watershed algorithm were explored. The Python-based computer vision library OpenCV [44] was used for these tasks.

The initially selected corner- and edge-detection algorithms were intended to recognize the boundaries of the blocks through the best-selected algorithm in order to be able to subsequently estimate, from this information, their total area (ABP) within each analyzed photo.

The first method tested was the classic and commonly used Sobel–Feldman operator, which is an edge-detection operator that uses convolutions with Sobel matrices to highlight gradient changes in an image [45]. The other method applied was the Canny edge detector, based on the computational approach developed by Canny [46]. This edge-detection algorithm uses a multi-stage approach, including Gaussian smoothing, gradient computation, non-maximum suppression and hysteretic thresholding, to achieve robust and accurate edge detection in an image. This method improved previous approaches, introducing the possibility of arbitrary edge detection, with good robustness to noise and, hence, minimizing errors associated with false edges.

Even though focused on corner detection, other approaches were also tested, like the classic Harris Corner Detector and the FAST (Features from Accelerated Segment Test) method [47,48]. The first one (the Harris Corner Detector) is a corner-detection method that evaluates the local variation in intensity in different directions using a covariance matrix to identify points of interest or corners in an image. The second one (FAST) is an image-keypoint-detection method that uses an accelerated segment test to quickly identify corners and points of interest, being especially computationally efficient.

The Oriented FAST and Rotated BRIEF (ORB), a method for detecting and describing keypoints in images that combines the efficiency of the FAST detector with the robustness of the BRIEF descriptor, was alternatively considered. This method was proved to be rotation-invariant and resistant to noise [49]. The ORB technique was used, in the current case, by applying the functions cv2.ORB create() and detectAndCompute() over previously grayscaled images.

The last method used for this initial approach consists in image segmentation by means of the watershed transform algorithm, included in the OpenCV library [44,50]. This watershed algorithm is based on the safe extraction of the background and foreground and then the use of markers to make the watershed run and detect the exact boundaries by using the function cv2.watershed(). As in the already mentioned corner- and edge-detection methods, the images require an initial grayscale transformation and subsequent binarization to make this method as accurate as possible.

Overall, it has to be stressed that the performance of the aforementioned methods depends on some features of the analyzed photos. Specifically, the ability of the edge- and corner-detection methods is a function of appropriate thresholding, scale, noise and inherent properties of the image, like contrast, color heterogeneity and lighting. Images with low contrast or with a predominant color may not be suitable, as observed for the set of images analyzed in this work.

2.6. The Segment Anything Method Algorithm (SAM) and Its Implementation

SAM is a state-of-the-art computer vision model developed by Meta AI [27] for object segmentation in images. The model is based on a deep neural network architecture that has been trained on a large dataset of labeled images. The training dataset includes over 11 million images from a variety of sources, such as the COCO, Pascal VOC and Cityscapes datasets. In addition to the large number of images, SAM has also been trained on a large number of segmentation masks. The training dataset includes over 1 billion segmentation masks, which have been carefully annotated by human experts. These masks provide the model with the ground-truth information needed to learn how to accurately segment objects in images.

From the three Vision Transformer (ViT) model levels [51] available in SAM (SAM-ViT-B (Base), SAM-ViT-L (Large) and SAM-ViT-H (Huge)), SAM- ViT-H was selected since it represents the most complete option and was established as ground truth for comparative purposes between the models, as presented by Zhang et al. [52].

For the implementation of the SAM algorithm, the size of the images was firstly reduced from the original size (1536 × 2048 px) to 768 × 1024 px, in order to improve the image processing time, reducing it to a third of the original one. It was checked that the size reduction did not affect the obtained results.

Once the selected model had been defined (SAM-ViT-H), we generated the segmentation masks for the selected image using the SamAutomaticMaskGenerator class from the Segment Anything library for Python. This class generates a grid of point prompts over the image and then filters low-quality and duplicate masks. The default settings for the arguments in this class were chosen for SAM with a ViT-H model. In our case, the parameter used for calibrating each image was pps (“points per side”). This parameter defines the number of points to be sampled along one side of the image, and the total number of points is pps².

After the segmented masks were obtained, representing the potential rock blocks detected in the picture, a function called show_anns was defined as Algorithm A1 (Appendix A). This function is a useful tool for visualizing the masks, as well as for calculating their total area in squared pixels (px²). The function takes as input a dictionary anns containing information about the masks, such as their segmentation masks, areas, bounding boxes and point coordinates.

The ABP (%) is calculated by comparing the total area covered by all masked areas (the sum of pixels within these areas) with the total squared pixel (px²) area of the image (Equation (1)).

$$ABP\left(\%\right)=\frac{total\_area}{image\_width\times image\_height}$$

(1)

To generate the visualization, the function first sorts the masks by their areas, from smallest to largest. Then, it assigns a color to each one according to a defined rainbow colormap scale from the Matplotlib library in the function map_to_rainbow, previously defined. As mentioned above, in addition to the visualization, the function also calculates the total area of the masks, which allows for the estimation of the percentage of blocks (ABP) in relation to the total image area.

The performance of SAM was evaluated through three metrics, namely, precision, recall and F1-score. These metrics (defined by Equations (2)–(4)) are based on the comparison of the segmentations performed by SAM and those considered “ideal” (in the present case, those corresponding to the masks manually created with the help of Label Studio software [53] for the analyzed photos). The comparisons relate three parameters: “true positives” (TPs), that is, positive examples (blocks) that were correctly classified as positive by SAM; “false positives” (FPs), negative examples (not blocks) that were incorrectly classified as positive by the model; and, lastly, “false negatives” (FNs), positive examples that were incorrectly classified as negative by the model.

The metrics were determined with the help of Python’s scikit-learn library [54], specifically, with the functions precision_score(), recall_score() and f1_score(), from the images previously transformed to grayscale and, finally, binarized.

$$Precision=\frac{TP}{TP+FP}$$

(2)

$$Recall=\frac{TP}{TP+FN}$$

(3)

$$F1-score=\frac{2\times Precision\times Recall}{Precision+Recall}$$

(4)

3. Results

This section provides a concise overview of the main results obtained from the application of edge- and corner-detection techniques as well as from the watershed algorithm. The results obtained from these denoted conventional methods are presented along with the outcomes recovered from the implementation of the SAM model in terms of the recognition of rock blocks. The proficiency of SAM in generating histograms for the distribution of relative block sizes is also addressed.

3.1. Results from Conventional Methods

Neither the results derived from the edge- and corner-detection methods nor those obtained from the application of the watershed algorithm were valid in terms of edge or corner recognition that would allow ABP calculations.

Figure 5 highlights the low performance yielded by some of the analyzed methods, namely, the Harris corner detector (Figure 5a), the ORB method (Figure 5b), the Canny edge detector (Figure 5c), the watershed algorithm (Figure 5d) and the FAST method (Figure 5e).

In the case of the Harris corner detector (Figure 5a), even though the algorithm could identify several points (corners) on the contrasted areas between dark blocks and the matrix, with this information it was not possible to estimate a value of the ABP. For the ORB method (Figure 5b), only some points of the image were recognized by the algorithm, and these were also insufficient for ABP calculations. Regarding the Canny edge-detector operator, it was not able to recognize any rock block but only some contrasted areas due to shadows, as shown in Figure 5c. A realistic ABP value could not be obtained through the application of the watershed algorithm either, as can be seen in Figure 5e, since only some shaded areas were recognized. Finally, the FAST method, although capable of segmenting the geomaterial, did not allow the determination of a reliable value for the ABP (Figure 5e).

3.2. Results from SAM

In this section, the set of original photos taken at the study area is presented in Figure 6, Figure 7 and Figure 8 together with the overlapped and isolated masks over a black background in order to show the block recognition ability of SAM.

Figure 6 presents photos corresponding to different outcrops from the Orellán Formation. As shown, the algorithm is able to segment a good number of blocks observable in the images, excluding vegetation and even the voids left by already released blocks (Figure 6b,c). Nevertheless, the rock hammer used as a scale is always segmented by the method, something that occurs also with other reference elements, like the panel shown at the top of Figure 7a from the Las Médulas Formation and the measuring rods in Figure 7c,d.

In Figure 7b–d, some photos corresponding to materials from the Santalla Formation are shown. These images present a higher contrast than the rest, due to the fact that they were taken in an underground cellar with particularly low lighting and high moisture conditions. This environment caused some darkening of several blocks. Despite these challenging conditions, the algorithm successfully detected a significant number of blocks. This is demonstrated by the masks, which are clearly displayed on the right-hand side of Figure 7b–d.

In the case of Figure 8, three photos from the Orellán and Las Médulas Formations showing apparently different levels of rock block-to-matrix ratios were selected. Although the rock hammer was also detected, the increasing number or presence of blocks is reflected in the results given by the algorithm, both in terms of the masks assigned and, consequently, through the estimated ABP values.

For each picture in Figure 6, Figure 7 and Figure 8, a value of the ABP was calculated by means of Equation (1) and is provided in

Table 1. Results given by the SAM algorithm in terms of ABP for each of the analyzed photos, referred to the corresponding position in Figure 6, Figure 7 and Figure 8.

Image	Location	ABP (%)
IMG_6648	Figure 6a	23.28
IMG_6662	Figure 6b	25.35
IMG_6666	Figure 6c	15.03
IMG_6673	Figure 7a	35.26
IMG_6701	Figure 7b	63.82
IMG_6710	Figure 7c	57.90
IMG_6716	Figure 7d	64.60
IMG_6720	Figure 8a	17.23
IMG_6750	Figure 8b	31.93
IMG_6762	Figure 8c	88.71

. The algorithm also allowed the analysis of size distributions of the set of detected rock blocks, whose size is expressed in squared pixels (px²). Following this line, a histogram representative of each photo could be obtained after the segmentation analyses (see an example in Figure 9 and the rest—Figure A1 and Figure A2—in the Appendix B).

4. Discussion

In this study, we assessed the effectiveness of various edge- and corner-detection methods and a watershed algorithm alongside an artificial intelligence-based image segmentation tool, SAM, in estimating the areal block proportion (ABP) from photos taken in diverse environments (surface and underground settings). Additionally, the study sought to provide estimates of the relative size distribution of rock blocks at various outcrops. This information is extremely useful for geopractitioners, both for choosing the most appropriate tunneling excavation method [13,55] as well as for performing reliable numerical analyses. In fact, quantitative data on the granulometric distribution can be used to more rigorously set up the geometrical model of the engineering problem at hand [7,10].

It should be taken into account that, despite the relatively limited number of analyzed photos, these are representative of different scenarios within the study area, encompassing both surface and underground scenarios and various lighting and contrast conditions, as well as diverse rock block contents, sizes, colors and distributions.

Conventional edge- and corner-detection methods have been found to be inappropriate for resolving the main objective of this work. These methods require a dedicated and time-consuming pre-processing task for each image, which often leads to questionable performance. Such an approach is particularly impractical in varied environments, such as those found in our photographs, where extreme color homogeneity and poor lighting conditions are common. It is worth noting that the image capturing process was carried out under slightly stringent conditions—the absence of a tripod, the use of a non-professional camera and no specific lighting control—in an attempt to represent a survey conducted by a practitioner with minimal experience, rather than in a professional manner, so as to analyze the performance of SAM in such situations.

Also relevant is the fact that in almost all the images a rock hammer/measuring rod appears, which was initially placed on the ground to gain a clear idea of the dimensions of the elements when the initial approach was based on conventional segmentation methods. These elements were also segmented by the SAM model, so that their areas were calculated. In any case, they represent, at most, 5% of the total sum, a fact that could become particularly relevant for those scenarios with low ABP values. Moreover, it is worth noticing that most of the rock blocks that were not detected were small enough (i.e., they were smaller than $5\%\sqrt{A}$, where A represents the square root of the total image area) not to influence the rock mass properties at the scales of engineering interest and can be considered, according to Medley [15], not to be “geotechnically significant”.

From the metrics recovered after the application of Equations (2)–(4) (

Table 2. Obtained values from the analyzed metrics for each image.

Image	Precision	Recall	F1-Score
IMG_6648	0.46	0.89	0.61
IMG_6662	0.61	0.69	0.65
IMG_6666	0.50	0.64	0.56
IMG_6673	0.66	0.86	0.74
IMG_6701	0.58	0.96	0.72
IMG_6710	0.70	0.88	0.78
IMG_6716	0.59	0.94	0.73
IMG_6720	0.54	0.56	0.55
IMG_6750	0.59	0.86	0.70
IMG_6762	0.82	0.96	0.88

) and at a first glance, it can be observed that the results corresponding to the recall parameter are systematically larger than those corresponding to precision. This fact means that SAM tended to erroneously identify elements in the photos as blocks (i.e., selection of the rock hammer and measuring rod). In terms of the recall parameter, the majority of the values retrieved from Equation (3) are higher than 0.80, except those representative of IMG_6662, IMG_6666 and IMG_6720. In these cases, this result may be explained by the color homogeneity presented by these images. From a general assessment, the F1-score values are generally higher than 0.60, with the exception of the already mentioned images.

In an attempt to compare the performance of SAM regarding ABP estimation with that carried out by visual inspection of images by experts, a survey was conducted with 32 people. They were required to indicate, based on visual assessment solely and from the photographs presented in this study, the percentage of rock blocks (ABP) observed. The survey results are presented as boxplots in Figure 10, together with the corresponding average results and the ABP values given by the SAM algorithm.

Figure 10 evidences that all the answers provided are more or less normally distributed, with a dispersion (i.e., interquartile range) that is rather significant in almost all cases and generally greater than 20% VBP. A particular variability in the expert answers was obtained for the photos showing more contrast, like IMG_6662, IMG_6666 and IMG_6750 (Figure 6b,c and Figure 8b, respectively).

Most of the representative results (especially for seven out of ten) from each set of responses (averages and medians) are similar to the outcomes determined by the algorithm, indicating the ability of this tool to produce results comparable to those obtained through a number of visual expert assessments.

5. Application of the Approach to Understand the Mining Method at Las Médulas

The mining method used at Las Médulas, as described by “Pliny the Elder”, was named Ruina Montium [37], and the mentioned author slightly cited the facts of excavating drifts, cutting the supports and ultimately the collapse or caving of the mountain, so it is still not clearly understood. Some initial interpretations suggested possible sudden flooding of drifts and the generation of a “water hammer” effect (hydraulic shock) to produce the collapse of the mountains [33], but recent studies indicate that the water pressures achieved when using the available technology in Roman times would not be sufficient to produce this effect [56].

The presence of vertical and sub-vertical cliffs and mining drifts parallel to them suggested the possible use of those galleries in combination with controlled or propped slope undercutting as a feasible approach, according to the Roman practices available at the time, to generate the Ruina Montium phenomena [56]. This type of mechanism, able to generate secondary toppling due to the natural undercutting of slopes (associated with erosion) and ultimately producing sub-vertical cliffs, has also been observed in natural slopes [57,58].

The ABP estimates obtained through the already described methodology were crucial for the mechanical characterization of the block-in-matrix materials, a necessary step for correctly implementing the mechanism of the mining method used by the Romans at Las Médulas in a 2D FEM model.

5.1. Characterization of the Matrix

Due to the complex nature of these materials, it was impossible to obtain representative samples of the bimunit (conglomerate) to test. However, in the Las Médulas Formation, some blocks from bands of finer materials that make up the matrix were collected, as they fell from cliffs due to erosive processes (Figure 11a). This allowed the preparation of small specimens that were subjected to various laboratory tests, including X-ray diffraction, particle-size distribution from sieve analyses, density, uniaxial compressive and tensile strength tests, as well as point load tests.

These tests indicated that the mineralogical composition of this matrix presents roughly 55% quartz, 35% muscovite, 7% kaolinite, and some iron and titanium oxides. In terms of soil classification, the matrix corresponds to a silty sand 90% of which is composed of particles in the range of 0.1 to 2 mm. The average in situ density was 1.82 g/cm². The uniaxial compressive (UCS_matrix) and tensile strength (DTS_matrix) of the matrix were averaged at 1.5 MPa and 0.09 MPa, respectively (Figure 11b,c), from tests on prismatic samples (the only possible geometry that could be obtained). The cohesion (c_matrix) was determined as 0.4 MPa from Equation (5) [59], considering a representative friction angle of the silty sand, ϕ_matrix = 34° [60]. It is also relevant to mention that when subjected to crumb tests, 3 cm edge cubic samples of the matrix material totally disintegrated in 10 min.

$$c_{matrix}=\frac{{UCS}_{matrix}}{2{\left(1+\sin {\varphi }_{matrix}/1-\sin {\varphi }_{matrix}\right)}^{0.5}}$$

(5)

Based on in situ observations of human-made heaps of boulders, locally named murias and which resulted from removing and putting them apart to concentrate the gold in the matrix, the maximum angle of repose (α) of these bulk blocks was observed to reach 48° (Figure 11d).

5.2. Characterization of the Bimrock Material

From observations of the area and by resorting to the aforementioned method, it was possible to estimate representative values of the ABP (%) for the coarser conglomerates at the Las Médulas Formation at various outcrops. These values are in the range of 60–65%, similar to the estimates previously obtained by visual observations made by experts.

Kalender et al. [4] define “A” as a parameter representative of the bond or contact strength between the clasts and the matrix. In the case of the studied materials at the Las Médulas site, the adhesion can be considered moderate and the clasts were observed to be sub-rounded but also with angular shapes sometimes. In this way, by following Kalender et al.’s [4] methodology, an approximate value of parameter A = 110 was considered (associated with the already estimated angle of repose for the blocks, α = 48°).

The information required to characterize the mechanical behavior—cohesion (c_bimrock) and friction angle (ϕ_bimrock)—of the geomaterials from Las Médulas was obtained based on Equations (6)–(8) by following the approach proposed by Kalender et al. [4]. For these calculations, the inputs estimated in Section 5.1 were used. Note that a conceptual, yet sufficiently representative for the present case, simplification [61] was made, assuming that ABP ≈ VBP (being VBP = 65%).

$${\varphi }_{bimrock}={\varphi }_{matrix}\left[1+\frac{1000\left[\left(\alpha /{\varphi }_{matrix}\right)-1\right]}{1000+5^{\left(\left(100-VBP\right)/15\right)}}\times \frac{VBP}{VBP+1}\right]$$

(6)

$${UCS}_{bimrock}=\frac{\left(A-A^{\left(VPB/100\right)}\right)}{A-1}{UCS}_{matrix}$$

(7)

$$c_{bimrock}=\frac{{UCS}_{bimrock}\left(1-\sin {\varphi }_{bimrock}\right)}{2\cos {\varphi }_{bimrock}}$$

(8)

Therefore, and according to Kalender et al. [4], the geomechanical parameters of the materials found at Las Médulas could approximately be estimated as $c_{bimrock}$ = 240 kPa and ${\varphi }_{bimrock}$ = 47.2°. The tensile strength of the matrix was estimated at 900 kPa through uniaxial tensile strength tests, and this value can be assumed to be representative of that of the rock mass.

The steepest slopes in the zone (Figure 2a), based on detailed laser-scanner surveying, dip vertically (about 90°) and are 40 m high. The highest one attains 80 m with an average slope of 70°. These slopes are currently stable, and the most prevalent natural instability mechanism could be circular failure, given their geomechanical nature. Therefore, the stability can be back-analyzed to compute the domains of cohesion–friction where these observed slopes should stand [62]. Considering the aforementioned assumptions and based on the circular failure charts proposed by Hoek and Bray [63], a graph with the corresponding cohesion–friction stable domains could be plotted (Figure 12), allowing the reader to appreciate how the estimated parameters are close to the limit of the instability domain but in the stable area, which ensures that the estimated parameters are consistent with field observations.

For these strength parameters and consistently with in-place observations, the generation of the slope collapse (Ruina Montium) will encompass, for a 30 m height slope, first the excavation of drifts parallel to it at two levels (5 and 20 m above the slope toe) and 10 m, on average, inwards. This would generate, as computed according to Kirsch’s equations [64,65] and as observed in place, vertical tension cracks growing from the drift centers upwards and downwards. Then, the slope would be undercut with a cut 2 m high and around 4 m long that should be timbered to ensure its stability. This having been done, and based on experience and timber behavior, these props would be removed, producing the collapse of the area. The hypothetical mechanism is illustrated in three steps in Figure 13.

A simple 2D FEM numerical model of the process described is illustrated in Figure 14, based on the software RS2 v. 11 [66], which shows the generated cracks and shear bands in terms of plastic shear strain, illustrating the potential separation of the large piece of conglomerate that would collapse after the indicated process.

The discretized area in the model encompasses an initial rectangle 200 m long and 80 m high. After excavation of a vertical 30 m high undercut slope, the presented results focus the center of the model. Dirichlet-type boundary conditions are provided, constraining horizontal displacements at both lateral sides and both vertical and horizontal displacements at the bottom. The code is run to attain an initial equilibrium before mining, representing an elastic stress state. Then, the vertical slope is excavated, the drifts are mined, the slope is undercut in two supported stages and, finally, the support is removed (Figure 14c).

Some kind of control based on experience is supposed to have been used in order to avoid hurt to the miners, even if the safety standards at the time were probably not too stringent.

6. Conclusions

The geomechanical behavior of block-in-matrix materials is strictly dependent on their volumetric block proportion (VBP). The quantity and size distribution of rock blocks can be estimated via 2D measurements, using stereological techniques, when sufficient data are available. ABP calculation is mostly performed manually by experts; therefore, it can be subject to human bias and error.

In this research, we attempted to use, on the one hand, several computer vision algorithms for the calculation of the ABP parameter. Conventional methods (the Sobel operator, the Canny edge detector, FAST, ORB and the watershed algorithm) showed slack performance due to the variety of the pictures analyzed and their dependence on numerous hyperparameters. On the other hand, the Segment Anything Model (SAM) achieved reliable detection of the different blocks in images from different sources, as analyzed from the determination of three metrics (precision, recall and F1-score). These results allowed the correct calculation of the corresponding ABP parameter in each case. This calculation was implemented as an ad hoc script, which is available in an open repository at https://github.com/iperezrey/PebDec (accessed on the 25 January 2024).

Additionally, these results were validated through comparison with the criteria of 32 experts who carried out a traditional visual analysis to estimate the ABP of each image. These advanced image analysis techniques have made it easier to obtain ABP estimates used to estimate geomechanical parameters needed to model excavation approaches used by the Romans 2000 years ago.

In conclusion, SAM emerged as the most fitting tool for ABP calculation, not only due to its general accuracy, but also due to its potential to reduce human bias as well as decrease the time spent by experts on repetitive tasks. The potential of the approach herein presented is thought to be especially helpful for applications encompassing complex geotechnical formations (such as bimunits) with photos of limited quality and diverse features.

In this particular study, the approach was applied for an estimate of conglomerate strength behavior that was used to research the practice of a mining method used by the Romans in the archaeo-mining site of Las Médulas in NW Spain. A mechanically sound hypothesis for Ruina Montium was in turn proposed, which needs confirmation based on the archaeological record.