Emilio Trigueros , Manuel Cánovas , Javier Arzúa and Manuel Alcaraz

Stability of an abandoned siderite mine: A case study in northern Spain

De Gruyter | Published online: March 23, 2021

Abstract

Bodovalle is an iron carbonate mine located in Gallarta (near Bilbao, northern Spain), which is currently in the closure stage. The deposit was first exploited as an open-pit mine and later as an underground mine. The underground mine currently has 40 large rooms with rib pillars, occupying an area of 2,000 m long by 600 m wide. Room depth is around 200 m. The main geotechnical incident occurred in the NW1N zone, where an over-exploited rib pillar in poor condition was partially removed, leaving a 70-m roof span that eventually collapsed in 1999, resulting in a surface crater measuring 60 m in diameter. The collapse was preceded by noises and movements detected outside the mine. The subsidence was improved by means of infilling and monitoring. In 2014, residents living in homes built over the abandoned mine rooms detected very loud noises similar to those of 1999. The article describes a stability review carried out using advanced numerical methods (finite element method and boundary element method) and the geomechanical parameters obtained from a back analysis of the 1999 collapse.

1 Introduction

Bodovalle is a semi-abandoned underground room-and-pillar mine located in Gallarta, a town with a population of approximately 10,000 inhabitants located in the Basque Country (northern Spain). The extracted ore was iron carbonate destined for the blast furnaces of Bilbao. The mine ended production in the 1990s and is currently in the closure stage.

From a geological perspective, the mineralization corresponds to Upper Aptian limestone and is formed of siderite – the most abundant ore – and ankerite to replace the limestone. The mineralized limestone is part of a formation consisting of calcarenites, limestone, marlaceous lime and sandy loams. The mineralized zone, occupying a graben with a sub-horizontal tabular shape, dipping from 10° to 30°, with a mean value on the central zone (where NW1N rooms are) of 15°, extends lengthways over 2000 m in the SE-NW direction and is 600 m wide on average. The iron concentration of the deposit is low (27%), but the company performed a process to increase its concentration to 40% before its commercialization.

Ore was initially extracted by means of open-pit mining and, once the pit reached the projected depth, an underground mine was opened using the room and pillar method. The underground mine, approximately 2,000 m long by 600 m wide, with large rooms and intermediate rib pillars, is formed of seven mining zones (Block 1, SE, NW Exp, NW Exp S, NW1N, NW1S and NW2; see Figure 1). Mining commenced in the NW Exp zone, with rooms 13 m wide, rib pillars 8 m wide and height 40 m, and very irregular due to the characteristics of the deposit. The over-burden exceeds 200 m, except in the area closest to the open pit, where the over-burden is 100 m.

Figure 1 General plan view of the Bodovalle mine, showing the faults conforming the graben and the layout of the rooms. The study area is the SE zone, marked in red, and the remaining zones are indicated in blue.

Figure 1

General plan view of the Bodovalle mine, showing the faults conforming the graben and the layout of the rooms. The study area is the SE zone, marked in red, and the remaining zones are indicated in blue.

The main geotechnical incident was a progressive pillar failure that occurred in the NW1N zone (Figure 1). The pillar, damaged due to a reduction in its width in order to increase the extraction ratio, finally collapsed in 1999 (Figure 2). The collapse, preceded by noises audible on the surface, resulted in a 60-m diameter crater on the surface, which was improved by backfilling (28% of granular backfill and 72% of cement-based backfill) and subsidence monitoring. A similar event occurred in 1986 in the same zone, when a progressive caving process reached the surface, resulting in the formation of a 35-m deep sinkhole, with a small-diameter opening that later grew into a crater measuring 30 m in diameter. The crater was improved through backfilling from the surface. The NW1N zone can, therefore, be considered to be almost completely backfilled and the only backfilled zone. In August 2014, local residents heard very loud noises similar to those of the 1999 collapse, and subsequent inspections found several rock blocks from the roof in room #8 in the SE zone.

Figure 2 Aerial photograph of the surface crater resulting from the 1999 collapse (image by the authors).

Figure 2

Aerial photograph of the surface crater resulting from the 1999 collapse (image by the authors).

This was the context in which the authors aimed to evaluate the subsidence risk associated with SE room #8, analyse surface land occupation and review the planned actions of the closure project. Used were numerical methods [finite element method (FEM) and boundary element method (BEM)] to analyse the stability of the SE zone, considering geomechanical parameters obtained from previous reports and fine-tuned by means of a back analysis of the 1999 collapse.

2 Methodology

Classical empirical pillar strength formulae are typically used as a first approximation to the design of pillars in coal mines [1,2,3,4,5,6], hard rock [7,8,9,10,11] and metal mines [12,13,14,15]. Empirical methods, many of which include the height-to-width ratio, are based on the uniaxial compressive strength of the rock, which, in combination with stress estimates derived from the tributary area method, can be used to calculate the factor of safety for pillars [16]. The rapid advance in computer tools and the mathematical development of FEM software and BEM software enables the determination of not only the possible pillar collapse but also the expected deformation type, stress distribution and the points where rock may fail due to either tension or shear stresses [17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]. Numerical methods have also been used to reduce the operational losses in room and pillar method [35,36] or to determine the original heterogeneous and anisotropic stress field [37].

Due to its origins and later deformations, the Bodovalle deposit is quite irregular and so the underground rooms are also irregular in shape and size. Empirical methods are based on regular pillars, however, and so consider that (i) the rock forming the pillars is homogeneous, (ii) the geometry of rooms and pillars is regular and (iii) there is constant vertical stress depending on the over-burden. However, these conditions are not fulfilled in the SE zone in Bodovalle: (i) room length varies, as it is adapted to the distance between local faults in most cases; (ii) the rooms are not usually in the same horizontal plane, which means that pillars have different face heights; (iii) the room height usually varies; (iv) some pillars combine different types of rock and (v) the over-burden is affected by the open pit.

For these reasons, this research was performed using advanced computational tools instead of empirical methodologies. Development of the models, i.e. determining the geometry of the excavations, the geological configuration of roofs and pillars and the geomechanical parameter values, was a complex task due to the inherent difficulties implied by the irregularity of the deposit and the semi-abandoned state of the mine.

Uniaxial compressive strength and rock quality designation (RQD) values were obtained from previous studies and the geological strength index (GSI) [38] was determined from RQD values. The GSI allows global rock mass strength to be calculated using an equation proposed by Hoek et al. [39]. The value of K (the ratio between horizontal and vertical in situ stresses prior to excavation) should fall in the interval between 0.8 and 2.5 according to the work carried out by Hoek and Brown [40] in metal mines around the world. Given the recommendation of Hoek and Brown and the relatively smooth orography, K = 1.0, indicating that horizontal and vertical in situ stresses were equal before the excavation works.

The study covered the following phases:

  1. Geotechnical characterization based on an analysis of previous studies, specifically, the geotechnical reports corresponding to the closure project. It is to mention that these reports performed laboratory tests on the different materials of the rock mass.

  2. Topographic survey of the excavations in the SE zone, consisting of 3D laser scanner and photographic reports.

  3. Field inspection of the open pit and of the NW1 and SE zones of the underground mine to compare geotechnical characteristics.

  4. Back analysis of the 1999 collapse to fine-tune the geotechnical parameters used in the FEM software (Phase2 V8.0 [41]).

  5. Stability estimates for the SE rooms, according to the conditions established in the closure project, using the FEM software and BEM software (Examine 3D, [42]).

  6. Evaluation of the subsidence risk for SE room #8, analysis of surface land occupation and a review of the planned closure project actions.

Since geometric parameters are keys to calculating room stability, the reliability of data from the previous works was checked. It was therefore necessary to topographically survey SE rooms #6, #7 and #8 to determine their geometry, relative positions and pillar positions.

The 3D laser scanner technology was selected as the most suitable approach, to the topographic survey, due to the difficulties of access and poor visibility and also for safety reasons (to minimize the time spent in the mine). An ultra-high-speed Leica Scanstation P20 scanner (available from the Mining and Civil Engineering Department of the Universidad Politécnica de Cartagena) was used, given its ability to scan up to 1 million points per second and its extended range of up to 120 m.

Some 60 million 3D points were scanned and the 3D point cloud was processed using Cyclone 9.0 software (developed by Leica). The digital geometry allowed the following conclusions to be drawn: (i) the widths of the three rooms and the pillars were practically identical to those documented in previous works; (ii) the extended range of the 3D laser scanner was insufficient to measure the total lengths of SE rooms #6 and #7 and (iii) the cross-sections of the rooms were not rectangular, the corners between walls and roofs were curved and the corners between wall and floor appeared funnel shaped (see Figure 9).

3 Geotechnical context

Field inspections pointed to good correlation between joint sets and large-scale faults. The orientation of the SE rooms was favourable to the stability of potentially failing wedges in the walls. Since neither the main faults nor the joint sets can open a way to the excavation, there was no risk of falling wedges defined by these discontinuities. Roof rock blocks formed from cayuela, a local name for a type of carbonated siltstone, had fallen since the mine was closed, attributed to production blasting that caused fracturing and mechanical decay of the first few meters of rock into the rock mass. This circumstance was taken into account through coefficient D in the Hoek–Brown failure criterion, which reflects the disturbance experienced by the rock mass as the result of blast damage and stress relaxation. D varies from 0 for undisturbed in situ rock masses to 1 for very disturbed rock masses [39].

As well as the assignment of individual mechanical strength values to the different rock types, the following were taken into account: (i) the impact of production blasting on roofs, (ii) changes in rock types, most specifically in pillars, (iii) different degrees of jointing in the pillars according to the analysed zone and (iv) the existence of mylonitized areas surrounding the main faults. For each type of material, we calculated the GSI [38] and, using RocLab [43], the generalized Hoek–Brown criterion [39]. This failure criterion yielded Young’s modulus (Em), average uniaxial compressive strength for the pillars (σcp) and the maximum tensile strength supported by the roofs (σtp). Post-failure behaviour parameters (ar, mr, sr and dilation) were calculated according to Crowder and Bawden [44].

3.1 Geomechanical parameters

The following material types were identified: hard, weak and host cayuela; hard, slightly altered and mylonitized ore; footwall limestone and hanging wall sandy loams. Tables 1–4 summarize the geomechanical parameters for the different materials, based on laboratory tests carried out for previous research and field surveys. The values were calibrated through a back analysis of the 1999 collapse, as the best in situ test material available. Figure 3 shows images of the two types of cayuela.

Table 1

Geomechanical parameters of the hard cayuela

Hard cayuela
γ (MN/m3) 0.026
GSI 70
σci (MPa) 60 σcp (MPa) 18.09
σti (MPa) 5 σtp (MPa) 0.52
Ei (GPa) 39 Ep (GPa) 25
ν 0.3
mi (bibliography) 12
Hoek and Brown fit
Elastic Plastic
D = 0
mb 4.11 mr 1
s 0.0357 sr 0
a 0.501 a 0.501
Dilatancy 0.2

*The subscripts i, p and r stand for intact specimens, pillar scale and residual, respectively.

Table 2

Geomechanical parameters of intact and blast-damaged weak cayuela

Intact weak cayuela Blast-damaged weak cayuela
γ (MN/m3) 0.026 γ (MN/m3) 0.026
GSI 55 GSI 55
σci (MPa) 60 σcp (MPa) 14.12 σci (MPa) 60 σcp (MPa) 8.1
σti (MPa) 5 σtp (MPa) 0.14 σti (MPa) 5 σtp (MPa) 0.06
Ei (GPa) 39 Ep (GPa) 20 Ei (GPa) 39 Ep (GPa) 20
ν 0.3 ν 0.3
mi (bibliography) 15 mi (bibliography) 15
Hoek and Brown fit Hoek and Brown fit
Elastic Plastic Elastic Plastic
D = 0 D = 0.8
mb 3.007 mr 1 mb 1.03 mr 0.5
s 0.0067 sr 0 s 0.0011 sr 0
a 0.504 a 0.504 a 0.504 A 0.504
Dilatancy 0.2 Dilatancy 0.2

*The subscripts i, p and r stand for intact specimens, pillar scale and residual, respectively.

Table 3

Geomechanical parameters of the intact and mylonitized ore

Intact ore Mylonitized ore
γ (MN/m3) 0.033 γ (MN/m3) 0.033
GSI 70 GSI 40
σci (MPa) 80 σcp (MPa) 22.48 σci (MPa) 60 σcp (MPa) 8.38
σti (MPa) 8 σtp (MPa) 0.83 σti (MPa) 4 σtp (MPa) 0.06
Ei (GPa) 60 Ep (GPa) 30 Ei (GPa) 20 Ep (GPa) 5
ν 0.2 ν 0.2
mi (bibliography) 10 mi (bibliography) 10
Hoek and Brown fit Hoek and Brown fit
Elastic Plastic Elastic Plastic
D = 0 D = 0
mb 3.425 mr 1 mb 1.17 mr 0.5
s 0.0357 sr 0 s 0.0013 sr 0
a 0.501 A 0.501 a 0.511 a 0.511
Dilatancy 0.2 Dilatancy 0.2

* The subscripts i, p and r stand for intact specimens, pillar scale and residual, respectively.

Table 4

Geomechanical parameters of the footwall limestone and of the hanging wall sandy loams

Footwall limestone Hanging wall sandy loams
γ (MN/m3) 0.026 γ (MN/m3) 0.024
GSI 70 GSI 40
σci (MPa) 65 σcp (MPa) 18.26 σci (MPa) 15 σcp (MPa) 1.758
σti (MPa) 7 σtp (MPa) 0.68 σti (MPa) σtp (MPa) 0.023
Ei (GPa) 59 Ep (GPa) 29 Ei (GPa) Ep (GPa) 2
ν 0.25 ν 0.25
mi (bibliography) 10 mi (bibliography) 7
Hoek and Brown fit Hoek and Brown fit
Elastic Plastic Elastic Plastic
D = 0 D = 0
mb 3.42 mr 1 mb 0.821 mr 0.821
s 0.0357 sr 0 S 0.0013 sr 0.0013
a 0.501 a 0.501 A 0.5 a 0.5
Dilatancy 0.4 Dilatancy 0

* The subscripts i, p and r stand for intact specimens, pillar scale and residual, respectively.

Figure 3 (a) Hard cayuela in the roof of SE room #8. (b) Weak cayuela near the collapse location (images by the authors).

Figure 3

(a) Hard cayuela in the roof of SE room #8. (b) Weak cayuela near the collapse location (images by the authors).

3.2 Back analysis of the 1999 collapse

The main objective of this analysis was to fine-tune the geomechanical parameters of the different materials. A FEM of the NW1N zone (Figure 4) was performed according to the characteristics of the collapse as documented in the technical report of the failure:

  1. Five 25 m wide by 50 m high rooms with four 20 m wide rib pillars and an average over-burden of 220 m.

  2. One rib pillar (formed of mylonitized ore) in poor condition by the 1980s was partially removed, leaving a 70-m roof span that collapsed in 1999.

Figure 4 Initial numerical model of NW1N zone rooms where the 1999 collapse occurred.

Figure 4

Initial numerical model of NW1N zone rooms where the 1999 collapse occurred.

Figure 5 depicts the collapse sequence, which follows the different stages of the excavation. After excavation of room #1, yielding rock appears in the limestone floor, in the upper part of the left wall and in the cayuela roof because of the damage produced by production blasting, although dilation prevents rock blocks from falling. The left wall of room #1 does not yield, apparently because of the mylonitized ore a few meters to the right (inside the rock mass), which, rather than loading that left wall, liberates energy in the form of displacement (Figure 5a). When room #2 is excavated, developments are similar. Again, deformation tends to concentrate in the mylonitized ore in the future pillar #2. The yielding occurs mainly in the wall of this future pillar #2 due to its poorer quality (resulting from mylonitized ore, blasting damage and over-excavation), which bears most of the deformation in order to accommodate stresses (the upper right corner of room #2 is vertically displaced by 1 cm, while a similar displacement of 0.8 cm occurs in the upper left corner) (Figure 5b). After excavation of room #3, pillar #2 deteriorates further and deformation continues (the upper right corner of room #2 is vertically displaced by 4.5 cm) to the point that yielding begins to occur in the previously healthy pillar #1 (Figure 5c). After excavation of room #4, further tension yielding points appear in the roof of rooms #2 and #3, pillar #2 is further vertically displaced (6.85 cm in the upper right corner of room #2) and yielding now extends to pillar #3 (Figure 5d). After excavation of room #5, pillar #2 is vertically displaced by 10 cm and the deformation leads to progressive damage of the pillar (Figure 5e). The pillar finally collapses, leaving a 70-m roof span, where yielding progresses upwards to reach the full height of the weak cayuela (70 m) (Figure 5f). The progressive collapse finally reaches the surface.

Figure 5 (a–f) Excavation and collapse sequence of NW1N zone rooms, explanation in the text.

Figure 5

(a–f) Excavation and collapse sequence of NW1N zone rooms, explanation in the text.

Since the FEM results were confirmed to faithfully reflect the technical report of the failure, the material properties, as based on previous geotechnical studies and the bibliography and fine-tuned by means of the FEM, were considered suitable to calculate SE room and pillar stability.

4 SE rooms stability analysis

Two different numerical methods were used to evaluate the stability of SE zone rooms and pillars: 2D FEM for different cross-sections, followed by 3D BEM of the entire zone.

4.1 2D FEM results

The FEM (Phase2 v8.0 [41]) was performed to explore the influence on stability of the excavations of different rooms. Note that this study only concerns itself with the final excavation, as the objective of the study was to evaluate the future risk of subsidence over the excavated rooms.

Since the geometry and locations of the rooms were different, four different cross sections of the SE zone were considered, as depicted in Figure 6, which shows the risky situation of the area located over rooms #3, #4, #5, #6 and #7 if any of these rooms collapsed.

Figure 6 Location and geometry of the SE zone rooms.

Figure 6

Location and geometry of the SE zone rooms.

In what follows, the strength factor (SF) – also known as the stress safety factor – which is calculated by dividing the rock strength (obtained from the used failure criterion) by the induced stress at every point in the mesh [41], while the strength reduction factor (SRF) obtained by means of the shear strength reduction technique; in this method, the strength parameters of a model are reduced by a certain factor (SRF), when the FEM model does not converge (it becomes unstable) the SRF is called critical SRF, being considered another way to calculate the safety factor of the model [41].

4.1.1 Northern cross section

The northern cross section goes through rooms #3 to #6 and through the northern end of rooms #2 and #7 (Figure 6). This section was selected for analysis because of the road and the children’s playground located over rooms #6 and #7. Figure 7 shows the FEM geometry and the materials. Note the different geometries of the rooms, going deep into the footwall of limestone in some cases and remaining far from the hard cayuela of the roof in other cases. This situation and the geometry are probably due to a lack of planning. Noteworthy is the extraordinary height of room #6 (84 m) and also the fact that, except for rooms #6 and #7, the rooms have ore in the roof.

Figure 7 Geometry and materials of the northern cross section of the SE zone rooms.

Figure 7

Geometry and materials of the northern cross section of the SE zone rooms.

Analysis of the FEM results (Figure 8) reveals the following:

  • The largest displacements (around 1.3 cm) are in the walls of room #6, while surface displacements are around 0.4 cm and occur approximately over room #3 (Figure 8a).

  • The SF is high in most cases, with the lowest values occurring in the pillars but do not develop completely between rooms. The lowest SF values (1.26) are below room #5 and in the left wall of room #7, with both tending to the lower part of room #6, while, following this trend, only a few elements show yielding in the floor of room #5 (Figure 8b).

  • High concentrations of principal stress (mostly vertical) occur in the left walls of the pillars between rooms #4 and #5 and between rooms #5 and #6, reaching values of 15 MPa (Figure 8c). The floor of room #6 also achieves this principal stress value (mostly horizontal), with the square corners of room #6 as the probable cause.

  • In applying the SRF, Phase2 points to a critical value of 1.36 (the strength parameters should reduce a 36% to get the non-convergence of the model, i.e. the failure of the rooms), with the failure surface occurring between the lower parts of rooms #6 and #7 (Figure 8d). According to the results of the model, if the SRF value reaches 1.4 (the strength parameters reduce a 40%), the failure could become generalized and affect the rest of the rooms, and if it is higher than 1.7 (the strength parameters reduce a 70%), the collapse could affect the surface.

Figure 8 Finite element method results for the northern cross-section of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

Figure 8

Finite element method results for the northern cross-section of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

4.1.2 Central cross section

The central cross section goes through rooms #2 to #8. This section was selected because the slope of the open pit may affect rooms #7 and #8, as will be explained below. Figure 9 shows the FEM geometry and the materials. Contrasting with the northern cross section is the similar geometry of these rooms, with rounded corners in the roofs and funnel-shaped floors. A 5-m (determined as a function of the tension strength of the intact rock and the peak particle velocity originated by production blasting [45]) layer of damaged cayuela was added to the roof of each room to take into account the damage produced by production blasting.

Figure 9 Geometry and materials of the central cross section of SE zone rooms.

Figure 9

Geometry and materials of the central cross section of SE zone rooms.

Analysis of the FEM results (Figure 10) reveals the following:

  • The displacements are small, reaching a maximum of around 1.4 cm in the roof of room #5. The shape of the displacements is very influenced by the open pit, which causes asymmetry of the loads and, therefore, of the displacements. Surface subsidence in the order of 0.8 cm occurs approximately over room #4 (Figure 10a).

  • SF reaches a value as low as 1.02 in the pillar between rooms #5 and #6 (Figure 10b). This low value occurs because there is no longer damaged cayuela in the roof. The asymmetric load (because of the open pit) combined with the great height of the pillars produces some buckling in the pillars, which ultimately causes yielding in some elements (Figure 10d).

  • Also because of the asymmetric load, the largest compression values do not occur in the pillars but in the ore base of the pillars between rooms #3 and #4 and between rooms #4 and #5 (Figure 10c). Similarly, the maximum tension values do not occur in the roofs, but in the pillars, reaching a value of −0.45 MPa in the pillar between rooms #6 and #7.

  • In applying the SRF, Phase2 points to a critical value of 1.24, with the failure surface occurring in the middle of the pillars between rooms #4, #5, #6 and #7. Note that rooms #5 and #6 will be backfilled according to the closure project, so this will ultimately improve the stability of these rooms.

  • Note also that, according to FEM results, room #8 does not play a key role in the stability of the cross section.

Figure 10 Finite element method results for the central cross section of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

Figure 10

Finite element method results for the central cross section of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

4.1.3 Southern cross section I

Southern cross section I goes through rooms #2 to #7 and correspond to the southernmost end of rooms #2, #3, #4 and #7 (Figure 6). In this cross section, the road is above room #2 and the remaining rooms are at the current limits of the open pit. Figure 11 shows the FEM geometry and materials. Also located in this zone is a unit of sandstone below some rooms and mineralized rock that seems to go through the limestone unit.

Figure 11 Geometry and materials of southern cross section I of the SE zone rooms.

Figure 11

Geometry and materials of southern cross section I of the SE zone rooms.

Analysis of the FEM results (Figure 12) reveals the following:

  • The displacements are small, at a maximum of around 1.1 cm in the roof of room #4 and 1 cm in the roof of room #5. The shape of the displacements is, again, very influenced by the open pit, which causes asymmetry of the loads and, therefore, of the displacements. Surface subsidence in the order of 0.7 cm occurs approximately over the pillar between rooms #3 and #4 (Figure 12a).

  • The SF reaches values of around 1.11–1.12 in the pillars around room #5 (Figure 12b). The lowest value in this section occurs in the contact between hard cayuela and limestone. Note that the closure project plans for backfilling of rooms #5 and #6, which will improve the stability of these rooms.

  • As in the central cross section, the largest compression values do not occur in the pillars but in the ore base of the pillars between rooms #3, #4 and #5 (Figure 12c), reaching values of 14.5 MPa. Similarly, maximum tension values occur in the pillar between rooms #5 and #6. The reason is, again, the asymmetrical loads caused by the open pit.

  • In applying the SRF, Phase2 points to a critical value of 1.25, with the failure surface occurring in the pillars between rooms #3, #4 and #5 (Figure 12d).

Figure 12 Finite element method results for southern cross section I of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

Figure 12

Finite element method results for southern cross section I of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

4.1.4 Southern cross section II

Southern cross section II comprises the southernmost end of rooms #5 and #6 (Figure 6). This cross section was considered because, since it is the nearest point to the slope of the open pit, there was some concern about how the rooms affected the slope and vice versa. Figure 13 shows the FEM model geometry and materials. No limestone appears in this cross section.

Figure 13 Geometry and materials of southern cross section II of the SE zone rooms.

Figure 13

Geometry and materials of southern cross section II of the SE zone rooms.

Analysis of the FEM results (Figure 14) reveals the following:

  • The displacements are small, at a maximum of around 0.7 cm in the roof of room #5 and 0.6 cm in the upper left corner of room #6. Deformation of the slope of the open pit becomes evident but is only 0.4 cm at the nearest point to room #6 (Figure 14a).

  • The SF reaches a minimum value of 1.22 in the upper left corner of room #6, where there is contact between hard cayuela and ore (Figure 14b).

  • The maximum compression stress occurs in the left base corner of room #5, reaching a value of 13 MPa (Figure 14c). Tension stress occurs in the remaining rock mass between the slope and room #6, reaching a value of −0.35 MPa.

  • In applying the SRF, Phase2 points to a critical value of 1.79, with the failure surface vertically connecting the roof of room #6 with the slope of the open pit (Figure 14d).

Figure 14 Finite element method results for southern cross section II of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

Figure 14

Finite element method results for southern cross section II of the SE zone rooms. (a) Total displacement. (b) Strength factor. (c) Major principal stress. (d) Strength reduction factor.

4.1.5 Effect of backfilling

As mentioned, the closure project indicates that rooms #5 and #6 will be backfilled. The same analyses as above were performed for fully backfilled rooms #5 and #6. The backfill was considered to be cohesionless granular material (mining rock waste) with a density of 2 t/m3 and a friction angle of 40°. Figure 15 and Table 5 depict the SRF results for the four cross sections, indicating great improvement after backfilling.

Figure 15 Strength reduction factor results after backfilling. (a) Northern cross section. (b) Central cross section. (c) Southern cross section I. (d) Southern cross section II.

Figure 15

Strength reduction factor results after backfilling. (a) Northern cross section. (b) Central cross section. (c) Southern cross section I. (d) Southern cross section II.

Table 5

Strength reduction factor before and after backfilling of rooms #5 and #6

Strength reduction factor Northern cross section Central cross section Southern cross section I Southern cross section II
Before backfilling 1.36 1.24 1.25 1.79
After backfilling 1.66 1.55 1.41 2.41

4.2 3D BEM results

Geomechanical 3D analysis allows the mechanical behaviour of the rooms as a whole to be considered, including the galleries connecting the different rooms. The objective is to identify zones with stress concentrations or excessive displacements that may act as a trigger for greater instability.

The software used was Examine 3D [42], a BEM that enables complex geometries to be analysed. In our case, some simplifications were necessary. Thus, in situ stresses were constants and K = 1, and only one material was considered for roofs and pillars, although the model was run for all the existing materials in the rock mass. The resulting BEM is shown in Figure 16.

Figure 16 3D boundary element method model of the SE zone and its rooms.

Figure 16

3D boundary element method model of the SE zone and its rooms.

4.2.1 Strength factor results

The Examine 3D software presents results graphically, so it is possible to visualize a perspective view with coloured SF contours. In this form of output, instability is marked by red and orange, while yellow, green and blue indicate stability.

In Figure 17, which depicts the 3D BEM results for the different materials, it can be observed that some yielding seems to occur in the galleries connecting the rooms, in limestone and in hard cayuela (Figure 17c and d, respectively) but decreases as the geomechanical quality of the material improves (ore and hard ore, Figure 17c and d, respectively). The yielding zones are small, nonetheless, so if a failure occurs, the surviving pillars will support the remaining rock mass. A horizontal yielding zone also occurs in the east wall of room #4 that, again, is small.

Figure 17 3D boundary element method strength factor results. (a) Limestone. (b) Hard cayuela. (c) Ore. (d) Hard ore.

Figure 17

3D boundary element method strength factor results. (a) Limestone. (b) Hard cayuela. (c) Ore. (d) Hard ore.

From the analysis, it can be concluded that no geomechanical problems in terms of stress are likely to compromise the general stability of the excavation.

4.2.2 Displacement analysis

The Examine 3D software can also estimate displacements of the boundaries of the model. Figure 18 shows the displacement for the different materials. Since the different materials have different elasticities, colour coding is not constant in all the graphics in Figure 18. Nevertheless, the computed displacements do not exceed 2 cm, corroborating the obtained FEM results.

Figure 18 3D boundary element method displacement results. (a) Limestone. (b) Hard cayuela. (c) Ore. (d) Hard ore.

Figure 18

3D boundary element method displacement results. (a) Limestone. (b) Hard cayuela. (c) Ore. (d) Hard ore.

4.3 Previous empirical results

The closure project, developed in 2004, includes stability analyses of the SE zone in order to propose reinforcement measures if required. The stability analyses were conducted using the Panek method [13], according to which rock strength is modified depending on pillar height, width and length in order to estimate the pillar strength. Comparing pillar strength and tributary load allows a factor of safety to be calculated for each pillar.

All SE zone pillars obtained factors of safety above 3 (mean 5.2). For the collapsed zone once backfilled, a mean safety factor of 3.2 was obtained, for a minimum value of 2.9.

Note that the empirical methods seem to over-estimate factors of safety in comparison with the factors of safety obtained using the numerical methods described above. The explanation for this discrepancy seems to be non-fulfilment of the requirement to apply pillar strength empirical methods.

5 Conclusions

This stability study of the Bodovalle deposit was performed using data from previous reports, the closure project, in situ geological observations, topographical surveys and back analysis of a previous collapse, which allowed the parameters used in the numerical models to be fine-tuned.

Because of the complexity of room geometry, the differing materials forming roofs and pillars and the presence of a nearby open pit, the current in situ stress is complex and causes asymmetry of the loads on the pillars. For these reasons, it was considered more appropriate to use numerical methods rather than empirical methods to estimate the loads on the pillars. The empirical methods over-estimated safety by establishing the factor of safety values that were more than double those calculated by the numerical methods.

The maximum stresses computed by FEM and BEM are around 15 MPa and −1 MPa for compression and tension, respectively. The combination of these stresses produces some very localized yielding zones, with safety factors of around 1.1. Nevertheless, the load that these yielding elements are unable to handle is supported by the surrounding elements. The BEM analyses do not reveal any important geomechanical problems that might lead to global failure. The projected backfilling of rooms #5 and #6 will improve the stability of the excavations. In any case, the safety levels that can be deduced from the low SF values should dissuade the implementation of infrastructures or activities that could put people at risk.

Towards the end of a mine’s life, all its geomechanical circumstances (instabilities, rock types, weathering grades, joint characteristics, etc.) are well known, making it possible to develop detailed numerical models for inclusion in a closure study that will guarantee the safety of the mine. Also meriting attention is the importance of taking into account the local deterioration of the rock.

References

[1] Greenwald HP, Howarth HC, Hartmann I. Progress report: experiments on strength of small pillars of coal in the Pittsburgh bed. USBM, R. I; 1941. p. 3575. Search in Google Scholar

[2] Steart FA. Strength and stability of pillars in coal mines. J Chem Metall Min Soc South Afr. 1954;54:309–25. Search in Google Scholar

[3] Salamon MDG, Munro AH. A study of the strength of coal pillars. J South Afr Inst Min Metall. 1967;67:56–67. Search in Google Scholar

[4] Salamon MDG. Stability, instability and design of pillar workings. Int J Rock Mech Min Sci Geomech Abstr. 1970;7(6):613–31. Search in Google Scholar

[5] Hustrulid WA. A review of coal pillar strength formulas. Rock Mech. 1976;8(2):115–45. Search in Google Scholar

[6] Bieniawski ZT. A method revisited: coal pillar strength formula based on field investigations. Proceedings of the workshop on coal pillar mechanics and design. Vol. 9315. Pittsburgh, PA: US Department of the Interior, Bureau of Mines, IC; 1992. p. 158–65. Search in Google Scholar

[7] Hedley DGF, Grant F. Stope and pillar design for the Elliot Lake uranium mines. Can Inst Min Metall Bull. 1972;65:37–44. Search in Google Scholar

[8] Sjöberg J. Failure modes and pillar behaviour in the Zinkgruvan mine. Proc. of the 33th U.S. Rock Mechanics Symposium, Santa Fe. Rotterdam: A. A. Balkema Publ; 1992. p. 491–500. Search in Google Scholar

[9] Hardy P, Agapito JFT. Induced horizontal stress method of pillar design in oil shale. XV Oil Shale Symp. Colorado School of Mines. Golden, Colorado; 1982. Search in Google Scholar

[10] Sheorey PR, Loui JP, Singh KB, Singh SK. Ground subsidence observations and a modified influence function method for complete subsidence prediction. Int J Rock Mech Min Sci. 2000;37:801–18. Search in Google Scholar

[11] Alejano LR, Arzúa J, Castro-Filgueira U, Malan F. Strapping of pillars with cables to enhance pillar stability. J South Afr Inst Min Metall. 2017;117(6):527–40. Search in Google Scholar

[12] Lunder PJ, Pakalnis R. Determination of the strength of hard-rock mine pillars. Bull Can Inst Min Metall. 1997;90:51–9. Search in Google Scholar

[13] Panek LA. Estimating mine pillar strength from compression tests. Trans SME/AIME. 1981;268:1749–61. Discussion 272, 2005–2009. Search in Google Scholar

[14] Von Kimmelmann MR, Hyde B, Madgwick RJ. The use of computer applications at BCL limited in planning pillar extraction and design of mining layouts. Proc., ISRM Symposium, Design and Performance of Underground Excavations. London: British Geotechnical Society; 1984. p. 53–63. Search in Google Scholar

[15] Krauland N, Soder PE. Determining pillar strength from pillar failure observations. Eng Min J. 1987;8:34–40. Search in Google Scholar

[16] Brady BH, Brown ET. Rock mechanics: for underground mining. Springer Science & Business Media; 2006. Search in Google Scholar

[17] Hatzor YH, Talesncik M, Tsesarskym M. Continuous and discontinuous stability analysis of the bell-shaped caverns at Bet Guvrin, Israel. Int J Rock Mech Min Sci. 2002;39:867–86. Search in Google Scholar

[18] Bètournay MC. Abandoned metal mine stability risk evaluation. Risk Anal. 2009;29(10):1355–70. Search in Google Scholar

[19] Li LH, Yang ZF, Yue ZQ, Zhang LQ. Engineering geological characteristics, failure modes and protective measures of Longyou rock caverns of 2000 years old. Tunnel Undergr Space Technol. 2009;24:190–207. Search in Google Scholar

[20] Sunwoo C, Song WK, Ryu DW. A case study of subsidence over an abandoned underground limestone mine. Geosystem Eng. 2010;13(4):147–52. Search in Google Scholar

[21] Helm PR, Davie CT, Glendinning S. Numerical modelling of shallow abandoned mine working subsidence affecting transport infrastructure. Eng Geol. 2013;154:6–19. Search in Google Scholar

[22] Lafont G, Gunzburger Y, Mitri HS, Al Heib M, Didier C, Piguet J-P. Influence of mining excavation on energy redistribution and rockburst potential. 23 World mining congress, Aug 2013, Montreal, Canada; 2013. Search in Google Scholar

[23] Ahmed SS, Al Heib M, Gunzburger Y, Renaud V. Pillar burst assessment based on large-scale numerical modelling. Proc Eng. 2017;191:179–87. Search in Google Scholar

[24] Malli T, Yetkin ME, Özfirat MK, Kahraman B. Numerical analysis of underground space and pillar design in metalliferous mine. J Afr Earth Sci. 2017;134:365–72. Search in Google Scholar

[25] Vallejos JA, Delonca A, Perez E. Three-dimensional effect of stresses in open stope mine design. Int J Mining Reclam Environ. 2018;32(5):355–74. Search in Google Scholar

[26] Zhang C, Mitra R, Oh J, Canbulat I, Hebblewhite B. Numerical analysis on mining induced fracture development around river valleys. Int J Mining Reclam Environ. 2018;32(7):463–85. Search in Google Scholar

[27] Konicek P, Ptacek J, Waclawik P, Kajzar V. Long-termn Czech experiences with rockburst with applicability to today’s underground coal mines. Rock Mech Rock Eng. 2019;52:1447–58. Search in Google Scholar

[28] Georgieva T, descamps F, Gonze N, Vandycke S, Ajdanlijsky G, Tshibangu J-P. Stability assessment of a shallow abandoned chalk mine of Malogne (Belgium). Eur J Environ Civ Eng. 2020;1–15. 10.1080/19648189.2020.1762752. Search in Google Scholar

[29] Malinowska AA, Hejmanowski R. Building damage risk assessment on mining terrains in Poland with GIS application. Int J Rock Mech Min Sci. 2010;47(2):238–45. Search in Google Scholar

[30] Malinowska AA. A fuzzy inference-based approach for building damage risk assessment on mining terrains. Eng Struct. 2011;33(1):163–70. Search in Google Scholar

[31] Al Heib M, Duval C, Theoleyre F, Watelet JM, Gombert P. Analysis of the historical collapse of an abandoned chalk mine in 1961 in Clamart (Paris, France). Bull Eng Geol Environ. 2015;74:1001–18. Search in Google Scholar

[32] Salmi EF, Nazem M, Karakus M. The effect of rock mass gradual deterioration on the mechanism of post-mining subsidence over shallow abandoned coal mines. Int J Rock Mech Min Sci. 2017;91:59–71. Search in Google Scholar

[33] Salmi EF, Malinowska AA, Hejmanowski R. Investigating the post-mining subsidence and the long-term stability of old mining excavations: case of Cow Pasture limestone mine, West Midlands, UK. Bull Eng Geol Environ. 2019a;79(1):225–42. Search in Google Scholar

[34] Salmi EF, Karakus M, Nazem M. Assessing the effects of rock mass gradual deterioration on the long-term stability of abandoned mine workings and the mechanisms of post-mining subsidence – a case study of castle fields mine. Tunnel Undergr Space Technol. 2019b;88:169–85. Search in Google Scholar

[35] Skrzypkowski K, Korzeniowski W, Zagórski K, Zagórska A. Adjustment of the yielding system of mechanical rock bolts for room and pillar mining method in stratified rock mass. Energies. 2020;13(8):2082. Search in Google Scholar

[36] Skrzypkowski K. Decreasing mining losses for the room and pillar method by replacing the inter-room pillars by the construction of wooden cribs filled with waste rocks. Energies. 2020;13(14):3564. Search in Google Scholar

[37] Ahmed SS, Gunzburger Y, Renaud V, Al Heib M. Initialization of highly heterogeneous virgin stress fields within the numerical modelling of large-scale mines. Int J Rock Mech Min Sci. 2017;99:50–62. Search in Google Scholar

[38] Marinos VIII, Marinos P, Hoek E. The geological strength index: applications and limitations. Bull Eng Geol Environ. 2005;64(1):55–65. Search in Google Scholar

[39] Hoek E, Carranza-Torres C, Corkum B. Hoek–Brown failure criterion-2002 edition. Proc NARMS-Tac. 2002;1(1):267–73. Search in Google Scholar

[40] Hoek E, Brown ET. Underground excavations in rock. London: Institute of Mining and Metallurgy; 1980. Search in Google Scholar

[41] Rocscience. Phase2 v8.0 finite element analysis for excavations and slopes. Toronto, Ontario, Canada: Rocscience Inc; 2011. Search in Google Scholar

[42] Rocscience. Examine3D – a 3D computer-aided engineering analysis package for underground excavations in rock, program and users’s manual, version 4.0. Toronto, Ontario, Canada; 2007. Search in Google Scholar

[43] Rocscience. RocLab – rock mass strength analysis using the Hoek–Brown failure criterion, user’s guide, version 1.0. Toronto, Ontario, Canada: 2002. Search in Google Scholar

[44] Crowder JJ, Bawden WF. Review of post-peak parameters and behaviour of rock masses: current trends and research. Rocnews. 2004;fall:13. Search in Google Scholar

[45] Trigueros E, Cánovas M, Muñoz JM, Cospedal J. A methodology based on geomechanical and geophysical techniques to avoid ornamental stone damage caused by blast-induced ground vibrations. Int J Rock Mech Min Sci. 2017;100(93):196–200. Search in Google Scholar

Received: 2020-12-14
Revised: 2021-02-09
Accepted: 2021-03-05
Published Online: 2021-03-23

© 2021 Emilio Trigueros et al., published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.