Introduction

In raw minerals industry, bulk exploitation methods dwell on reducing mining costs by boosting production rates. For this purpose, gravity induced fragmentation and in-situ leaching techniques represent potential alternatives to the conventional drilling and blasting cycle. Although bulk mining methods provide operational and financial prospects, their viability is restricted by geochemical and geomechanical constraints leading to serious strata control issues. For example, the chaotic fracture network and irregular rock flow patterns render caving methods risky, while solution mining is restricted to chemically compatible minerals. In addition, extensive production openings with no roof support or backfill present formidable strata control challenges, including surface subsidence and sudden collapse. In this respect, an ideal mining method should balance safety and economic viability while concerning for the environmental and social implications (Nelson 2011). Statistics assert that increasing demand for raw minerals will rapidly deplete the shallow resources. For example, the total share of underground mining in Türkiye increased from 27% in 2006 (Eros and Candelario-Quintana 2006) to 31% by 2014 (Mining Intelligence 2014). In USA, surface mining decreased by 15% and underground mining increased by 14% between 2011 and 2015 (Osborne 2020). This trend highlights the escalating potential of strata challenges inherent to underground mining operations.

Advances in roof support systems and caving mechanics have improved the efficiency of longwall mining compared to the naturally supported open stoping techniques. However, caveability and groundwater flow increase the risks associated with sustainable production and occupational safety. Regarding these drawbacks, numerous studies have focus on roof control and pillar stability (Ünal 1983; Zipf 1999). Recent research has concentrated on the physical modeling (Ghabraie et al. 2015; Moradi et al. 2015) and computational simulations (Singh and Singh 2010; Chen et al. 2021) of longwall operations, exploring the effects of multiphysical factors (such as water and humidity) on caving and pillar design (Pechmann et al. 1995; Zipf and Swanson 1999; Liu and Liu 2021). Similar to the longwall practice, regular roof caving is critical in solution mining to sustain the global mine stability. However, deadweight of the roof strata may not always be sufficient to induce roof caving on stiff and thick geological units (Peng and Chiang 1984). Instead, stress concentration will increase the likelihood of sudden collapse. If the stiff layer takes place on the main roof, the weak bonding strength between the strata may allow the immediate roof to bend down into the opening but tensile stress accumulation will still pose a risk in the main roof in terms of deformations and seismic events (Stefanko 1983). Considering sedimentary deposits may be subject to various stratifications, long-term global stability requires concerning about the thickness and stiffness of geological layers in the roof.

Destressing is a stress management technique that mitigates mining induced stress accumulation. Destress blasting, destress drilling, water infusion, and hydrofracking are alternative techniques of inducing fractures within the stiff strata (Roux et al. 1957; Gu et al. 2019). Mining industry has long been taking advantage of blasting to create a stress shadow as a result of reduced stiffness in destress panels (Vennes and Mitri 2017). In deep hard-rock metal mines of South Africa, destress blasting has been performed on field scale (Tang and Mitri 2001; Konicek et al. 2011; Vennes et al. 2020). Also in longwall hard coal mines, several cases approve the viability of destress blasting to reduce burst potential in coal face (Konicek and Schreiber 2018; Wojtecki et al. 2020). Although efficiency is concerned due to relying on propagation of pre-existing fractures, field applications validate the method (Toper et al. 1997). Table 1 presents some example cases where destress blasting applied on the field for stress management (Vennes and Mitri 2017).

Table 1 A review of destress blasting cases (Vennes and Mitri 2017)

Aligning with stability assessment, computational modeling has a flexible framework for parametric analysis, which helps to determine the impact of destress blasting prior to its execution. However, the simulated stress and deformation responses are contingent on the numerical method employed. Comparing different numerical approaches Sainoki et al. (2016) demonstrated that traditional modeling techniques may yield overly optimistic results, whereas novel models simulate the damage more precisely. Despite the challenges of computational methods, the focus on static and dynamic simulations remains undiminished when designing and evaluating the performance of destress blasting (Momoh et al. 1996; Miao et al. 2022). Recently, advances in computational power and discontinuum based numerical codes have led to utilization of Discrete Fracture Networks (DFN) for improving model precision (Saharan 2010; He et al. 2017).

It is a common fact that solution mining decreases operational cost with bulk production. However, stiffness variations on the unsupported roof strata have potential to create uncertainty in terms of deformations. Similar to longwall mining, induced stress accumulation is a reason for sudden collapse. Common consequences are loss of production units and irregular hydraulic conductivity. While previous studies have explored destressing techniques in underground stoping methods and longwall mining, a clear understanding of destressing in solution mines remains elusive.

This research compares the stress relaxation of two generated destress blasting patterns using computational models. Firstly, ribbon pattern was simulated to determine whether stress relaxation can be forced in-between small-scale destressing panels by inducing stress concentration. Later, models of borehole pattern were investigated to test the potential of replacing large destress panels while maintaining effective stress relaxation. Parametric studies were conducted with a distinct element code and stochastic DFN models. Ribbon pattern was simulated for 30 m, 60 m and 90 m off-set distances and borehole patterns covered straight and staggered hole distributions. In advance, effect of destressed span was simulated for multiple stiff layers. Stress concentrations and the yielded region around a representative production unit were investigated by stress and deformation analyses. The study outcomes propose a scalable layout for destress applications under similar conditions.

Background of in-situ solution mining

In-situ solution mining diverges from the other leaching applications by dissolving the orebody in-place and recovering the solution via boreholes (Schlitt 2011). Although landmark is not directly disturbed due to the extraction of ore, this method is not classified as an underground mining technique (Boshkov and Wright 1973; Morrison 1976; Hartman 1987). Solution mining is employed to produce evaporites such as potash, nahcolite, and trona (Hardy et al. 2003; Nielsen et al. 2004) as well as rock salt and precious metals like copper, manganese, and uranium (Hore-Lacy 2004; Chamberlin 2009). Its history dates back to the early 1940s for trona mining (Pike 1943). Waste management is more effective in reducing environmental impacts compared to conventional mining since waste rock is left in place. Horizontal and vertical wells are used to reach and access the deposit. Hot solvent injection dissolves the minerals and forms the ‘brine’ that is later pumped out from the vertical wells. Various production layouts were developed regarding geology, geomechanical properties, production rate, and constraints like presence of an aquifer. Figure 1 illustrates the single-well, dual-well, or multiple-well layouts and variations of the multiple well method used in solution mining.

Fig. 1
figure 1

Single-well (a), dual-well (b), and multiple-well (d)(e)(f) layouts of in-situ solution mining

Field experience points out that chemical and mechanical limitations determine the production rate. For example, sodium bicarbonate blinding in trona is a processing challenge that arises when the bicarbonate ions generated during dissolution resettle on the exposed surface of the ore, impeding the extraction process. Impurities present between the thin and intercalating seams may also obstruct the extraction. Besides, the roof caving is a critical factor for sustainable production. As brine is extracted, the overlaying geological units are expected to yield and collapse into the unsupported production opening, called ‘cavern’. However, sedimentary rocks, formed by the accumulation of different minerals, exhibit variations not only in their mineralogical and physical properties but also in their geomechanical characteristics. While weak rocks can easily bend and fail, the stiff layers tend to bear more stress and pose a risk of sudden collapse. Case histories from room-and-pillar and in-situ solution mines point out that there is an association between the sudden collapse and strong overburden strata (Dowding and Andersson 1986; Ma et al. 2022). The strong roof can be helpful for creating large caverns; however, stress accumulation may disturb the rock mass on a larger span including the neighboring units (Daupley et al. 2005). Several occasions of serious roof collapses in solution mine caverns in the USA, France, Italy, and China have led to studies on cavern stability (Weishen et al. 2011; Xiao et al. 2016; Zhang et al. 2018, 2023; Liu et al. 2020). To mitigate the severe impacts of sudden collapse on both the production caverns and the surface topography, extensive research has been conducted on active and passive strata control methods.

Besides, challenges related to multi-cavern operations include pillar instability, management of aquifers and control of surface subsidence. Under these circumstances, continuous field monitoring is essential for an effective strata control. As solution mining provides no human access to the underground operation, the horizontal and vertical surface displacements are monitored by conventional field surveying methods or sophisticated satellite measurement techniques and used for calibration of geomechanical models (Morgan et al. 2018). In addition, micro-seismic monitoring is employed to trace post-mining effects of solution mining in the Netherlands and France (Kroon et al. 2003; Contrucci et al. 2011). Recent studies mostly focus on post-mining utilization of solution caverns for the storage of natural gas (Li et al. 2022), carbon dioxide (Shi and Durucan 2005), and hydrogen (Abreu et al. 2023).

Research methodology

This study employs numerical simulations to evaluate alternative destress blasting patterns in solution mining. Regarding the common geological and structural features involve multiple beds with contact relations, models were implemented using a Distinct Element Method (DEM) code. Stratification, geomechanical characteristics and cavern geometry were adopted from a real solution mine to create a hypothetical model for thick trona seams. Mechanical inputs were further calibrated based on the surface subsidence history of the benchmark mine. In order to investigate the effect of destressed span, stiff strata thickness was modified in hypothetical model. Once validating the ‘Discrete Fracture Networks (DFN)’ are capable of simulating the stress relief, three different ribbon patterns with 30 m, 60 m and 90 m off-set distance and two different borehole patterns in straight/staggered order were designed to investigate their performance. Each distinct production stage was imitated and the development of mechanical stress and deformations were modeled. Destressed region was represented with destress panels containing a fracture network in the ribbon pattern. For the borehole patterns, the fracture network was assigned around the close vicinity (the plastic region) of coarsely distributed boreholes. Both ribbon and borehole patterns were implemented separately in a thin and thick stiff layer in the overlying strata. In total, ten models were computed and performance of each method was evaluated based on stress and deformation analyses. Figure 2 presents the flow of the study.

Fig. 2
figure 2

Study flowchart

Computational simulation of a solution cavern

The large underground spaces produced by solution mining induce stress and deformations on the overlaying strata and around the opening. The mechanical response of the rock mass can be analyzed by physical models, field monitoring of a trial production, and geomechanical simulations. Although the laboratory scale physical models provide an overall understanding of stress flow, it is challenging to determine a rock-like material to represent the strata. Also, simplification of the structural and topological characteristics may dramatically affect the mechanical indicators. On the other hand, test production is a costly alternative as it requires field-scale operation and in-situ monitoring of stress and deformations. Under these circumstances, numerical modeling is the best alternative to observe various production sequences with parametric analyses. Once the rock mass behavior under various loading conditions is specified using laboratory or in-situ rock mechanics tests, computational simulations can be used to examine the strata movements in any of the production stages.

This section presents the numerical model of a solution cavern with stiff layers above the thick trona deposit. The model geometry, boundary conditions and the mechanical inputs are tested on a preliminary model to simulate the stress and deformations developing on the stiff layers. Various destressing patterns were implemented using the three-dimensional distinct element code, ITASCA 3DEC, to determine the most viable setting.

The model layout and inputs

A hypothetical solution mine model was generated based on the field and laboratory studies of a local production region in Beypazarı trona mine. The thickness of the trona seams and the mechanical inputs of two of the main roof strata were adjusted to create a model geometry with thick seams and stiff layers. Investigating the drill core samples, eight structural domains were determined and each of them were designated with a prefix ‘T’ and an incremental numerical tag starting with T1 on the surface and ending with T8 on the deepest layer. Table 2 shows the numerical model inputs.

Table 2 Rock mass Geomechanical parameters

‘Linear-elastic’ and ‘Mohr-Coulomb Strain-Softening’ models were employed to define the elastoplastic constitutive equations for the continuum blocks of evaporitic rocks. The material model was calibrated through laboratory uniaxial and triaxial loading tests on core specimen. Uniaxial monotonic loading tests under strain control yielded a peak stress to residual stress ratio of 3/2. In triaxial compression tests, only the peak strength was captured under confinement. However, the residual mechanical properties were determined regarding the experimental findings of Wang et al. (2018) on similar rock types. Accordingly, a linear trend between the residual normal stress and confinement was embraced with a maximum degradation of 10% from the peak stress versus confinement relationship. On the other hand, the contact surface mechanical properties of the natural surfaces were determined with laboratory direct shear tests and implemented in numerical models with the ‘Mohr-Coulomb sliding criterion’. For the destressed zone, normal and shear stiffness were reduced to 0.1% and the shear strength parameters were assumed to represent an idealized case of full separation and minimal interaction between blocks. Although assuming a full joint aperture with no contact in fracture surfaces may not fully account for the residual frictional resistance arising from asperities or dilatancy effects, this simplification was implemented to focus on the effect of destressing pattern. The contact mechanical properties are given in Table 3.

Table 3 Contact mechanical properties

In determining the model scale, a global-local modeling approach was embraced regarding the computational efficiency. Trial simulations with various model dimensions, grid densities and fracture networks (considered in terms of the number of distinct blocks and contact surfaces) have confirmed that a local model of greater detail supported with boundary conditions from the large scale is more suitable for the current study. Therefore, the initial model with a coarse grid was implemented with dimensions of 600 m in the x-axis, 1500 m in the y-axis and 500 m in the z-axis. The velocity boundary conditions were set on each face to completely fix the motion except the top plane imitating the topographical surface. The gravity loads were implemented with a hydrostatic stress condition that stands for around 8 MPa in-situ stress at the mining levels. In the following stage, stress boundary conditions (obtained from the large-scale model) were defined on a small-scale model with dimensions of 370 m in the x-axis, 800 m in the y-axis and 400 m in the z-axis. Also, the grid density was increased around the cavern. T7 layer presents the evaporite with a thickness of 60 m. The deposit was modeled to involve an upper and lower seam divided by a thick interburden. The seams and the interburden have an equal thickness of 20 m. The elevation of the lower seam is between 30 m and 50 m and the upper seam is between 70 m and 90 m. As the study aims to represent a common strata for solution mining, typical structural features were implemented in the upper and lower seams by discretizing the layers with vertical planes. In this manner, the whole trona deposit was divided into 15 sections in x-axis and 30 sections in y-axis. A total of 450 continuum blocks were created with approximately 24 m and 26 m in x and y axes. On the other hand, the roof strata and the interburden between the trona seams were not discretized due to the large model dimensions and increasing computational cost. Based on the case history, a rectangular cavern geometry was embraced with 70 m in width and 200 m in length (Fig. 3).

Fig. 3
figure 3

The model layout and a cross-section view showing a solution cavern and the simulated production stages

The major principal stress (σ1), the minor principal stress (σ3), and the vertical displacement were traced on the history points within the overlying strata for investigating the effect of destressing on T3 and T4 stiff layers. The history points were located right above the cavern, in each layer and along the centerline with 20 m spacing in between. The models were computed on a time-marched iterative manner and production was simulated in eight stages.

Analysis of mining induced stresses and strata displacements

In this study, computational models were used to evaluate the performance of various destressing strategies by comparing;

  • post-mining stress field without destressing,

  • post-mining stress field with destressing,

The initial motivation was to simulate the pre-mining stress field in a typical strata where solution mining is viable. A hypothetical strata model and production plan was modeled with real geomechanical characteristics obtained from Beypazari strata. Although in-situ stress ratio has influence on the induced stresses, it was fixed to hydrostatic field loads to focus solely on the effect of different destressing schemes. The model was calibrated to align with the expected pre-mining stress field based on the depth and equal horizontal and vertical stresses.

The operational experience of solution mining points out that the thin clayey interburden within the trona seams bend into the cavern and fill the void. To mimic this mechanism, the extraction steps were continued until the roof yields down and contacts with the floor. In the following stage, another trona layer was removed from the roof and this procedure repeated for each sub-layer. The simulation involves eight production stages in which the stress build up and vertical displacements were traced on the history points. The central history points in each layer above the cavern were selected along the cavern to compare the stresses.

Figure 4(a) shows the history locations along with the major principal stress contours after the first production stage. In Fig. 4(b), following the completion of the lower seam after the fourth production stage, an increase in the major principal stresses can be observed around the cavern. Figure 4(c) shows the stresses accumulation on T3 and T4 stiff layers after exploiting the upper trona seam. The major principal stress history of the overlaying strata shows that the failure of the T6 layer occurs immediately after the first production stage (Fig. 4(d)). After the third stage, the T5 layer fails. Also, the stress accumulates on T4 until it completely fails in the fourth production stage. The stress history proves that the stiff layers continue to bear load when the lower strata fail and bends into the empty space.

Fig. 4
figure 4

Evolution of the stress history in strata above the solution cavern

The vertical displacement histories of the overlying strata above a solution cavern are presented in Fig. 5, where T5 and T6 layers fail and settle drastically but no significant movement can be observed above T4. The maximum displacement in thin and relatively less stiff T4 layer is by 1.31 m. The displacements in the upper strata are 1.15 m in T3, 1.09 m in T2 and 0.98 m in T1.

Fig. 5
figure 5

The vertical displacement history throughout the production stages

The simulations prove that the stiff layers above the mining level do not move synchronously even in production of a single cavern. Uneven stress distribution may potentially threaten the global mine stability when multiple caverns are operated. Also, the stress history implies the potential of violent failure due to high stress concentration. So, destressing in the stiff layers is necessary to sustain the ground control.

Two destressing alternatives: ribbon and borehole pattern techniques

Considering economic and technical viability, two different destressing methods, which are ‘ribbon and ‘borehole’ patterns, were investigated. The ribbon pattern technique creates regularly spaced panels with an intense fracture network. Its objective is to produce weak regions to promote stress concentration in between the ribbons and induce destressing without any extra effort. In constructing the ribbons, hydraulic fracturing or blasting methods can be employed in closely spaced boreholes regarding the specific site conditions and the required fracture density. On the other hand, borehole pattern technique covers a stress accumulation zone with discrete boreholes arranged in a parallel or staggered layout and farther spacing compared to the ribbon technique. Although hydraulic fracturing and destress blasting methods also comply with the borehole technique, fracture density is likely to be compromised due to distance between boreholes. This study does not intend to simulate the fracture propagation but only investigates the effects of the induced fracture network on the rock mass. The study flowsheet is given in Fig. 6.

Fig. 6
figure 6

Study flowsheet

In the ‘ribbon’ pattern, the area above the cavern on the stiff layer was divided into multiple slices hosting a fracture network implemented by a stochastic DFN model. The method relaxes the whole stiff layer from bottom to top. The width of the ribbon was kept equal to the cavern width aligning with the cavern borders. Three different spacings (30 m, 60 m and 90 m) along cavern’s long axis were selected for evenly spaced ribbons within the stiff layers T3 and T4. While 4 ribbons fit into the roof strata with 30 m ribbon spacing, this number reduces down to 3 for 60 m spacing and 2 for 90 m spacing. Interval of the ribbon spacing was determined with common sense considering the cavern length, development time constraint and economic viability. Each ribbon is a panel of weakened rock mass with dimensions of 70mx10mx75m (52500m3) in T3 layer and 70mx10mx50m (35000m3) in T4 layer. In other words, 30 m ribbon spacing aims to relax 20% of the roof strata with stress accumulation, while this ratio reduces down to 15% for 60 m ribbon spacing and 10% for 90 m ribbon spacing. The models were coded by the prefix ‘RN’ designating the ribbon and the spacing separated by an underscore (RN_30, RN_60, RN_90) (Fig. 7). The models were coded with ‘RN’ standing for the ribbon and the spacing separated by an underscore. Comparing the stress history following the relaxation of the stiff strata, the best spacing was found.

Fig. 7
figure 7

The model body, dimensions, and pattern layout of RN_30 (a), RN_60 (b) and RN_90 (c)

The borehole pattern similarly applies on the same layers but the DFN relaxes only a small volume of rock. The boreholes were distributed with an equal spacing within the cavern dimensions and placed only above the cavern. Around each borehole a fractured region of 4 m in diameter was set to imitate the plastic region that can be implemented by destress blasting. The “straight” and “staggered” patterns were studied to check the effect of borehole pattern on the stress relaxation. In the staggered pattern (BH_STG), the row of boreholes was placed with 20 m spacing where the following row of the boreholes were placed in the middle of the former boreholes (Fig. 8 (a) & (b)) or simply in the staggered formation. In the straight pattern (BH_STR), the boreholes were placed symmetrically and distributed equally just over the cavern with 20 m spacing in x and y dimensions (Fig. 8 (c) & (d)). Both the ribbon and the borehole patterns were implemented individually on T3 and T4 layers to reveal the effect of destressing depth and thickness of the relaxed layer along with the patterns.

In both of the methods, a set of discrete, planar, and finite sized discontinuity surfaces were formed in the rock mass. Fracture density within the destressed zones was modeled using a stochastic DFN model, which simulates the fracture network based on input parameters such as mass-density (P32) of 0.25 and a spatial distribution. For the fracture size distribution, the power law was used as it is widely accepted for induced fractures. The fracture size limits were set between 50 m and 100 m. The position distribution of the fractures was defined to be uniform in 3D. For fracture orientation in terms of dip angle and dip directions, uniform distribution was selected.

Fig. 8
figure 8

The borehole pattern dimensions, distance, and the model body of the BH_STG (a) and BH_STR (b)

Results & discussion

Given the problem constraints are thick & stiff roof strata, and the operational & economic viability, the numerical models were designed with the following objectives:

  • determine the best spacing for the ribbon pattern in the upper stiff layer (T3) by comparing 30 m (RN_30), 60 m (RN_60) and 90 m (RN_90) distance between the destressed regions,

  • identify the best borehole pattern in the upper stiff layer (T3) layer by comparing a staggered (BH_STG) and straight (BH_STR) drillhole pattern,

  • compare the overall performance of ribbon and borehole destressing patterns in the upper stiff layer (T3),

  • analyze the effect of destressing in the upper (T3) and lower (T4) stiff layers by comparing the best practices of ribbon (RN) and borehole (BH) patterns,

  • assess the best destressing pattern relative to the model simulated without destressing (referred to as woDESS).

The major principal stress (σ1) history in Fig. 9 showed that both RN_30 and RN_60 patterns forced stress concentration between ribbons by around 10%. This is due to the fact that stress flow shifted from the weakened ribbons and focused on the undisturbed regions. Manipulation of in-situ stress field is another reason that provides stress relaxation. Fracture network diminishes the confining field stresses around the ribbons and that will result in degradation of the rock strength in between ribbons. Simulation shows that originally the maximum principal stress was around 6.5 MPa in the stiff stratum. However, the RN_30 pattern increases the concentration up to 7.1 MPa in between ribbons while it is 7.2 MPa for RN_60 pattern in the initial production stage. As production expands, the stress level drops gradually. This indicates weakening of the stratum without posing a sudden collapse risk. On the other hand, the stress level was observed to be around 6.2 MPa for the ribbon spacing of 90 m (RN_90). Under these circumstances, 90 m spacing considered insufficient for relaxation. Another remark was relaxation in the lower stiff layer (T4) as the upper stiff layer (T3) was destressed (Fig. 9). It was by 17% for RN_30 case, 21% for RN_60 case, and 23% for RN_90 case.

The simulations show that the stress accumulation between ribbons in the upper stiff layer (T3) leads to the relaxation of the whole stratum. Since the objective was to control the sudden collapse risk by weakening the stiff layers, the ribbon pattern was proven to be useful. With plastic simulations, the yielded elements and yield types were investigated and given for the final production stage (Fig. 9). The yield zone was observed to progress towards the topographical surface as the trona seams were produced. In addition, the effect of destressing in the stiff layers can be noticed. The models indicate a correlation between the fracture density and the ribbon spacing. Densely spaced ribbons increase deformations and stress relaxation.

In summary, the model interpretations point out that RN_30 and RN_60 are the most effective spacings for the ribbon pattern. There is no significant variation (around 2%) between them in terms of stress manipulation. However, RN_90 cannot provide an intense relaxation. Another conclusion is about the operational difficulties in practicing the induced fracture network. Considering the dimensions (70 m in width and 75 m in height), and the number of destress ribbons fitting on the roof of a single production cavern (four ribbons with RN_30, three ribbons with RN_60), the RN_60 pattern was deemed the most suitable pattern.

Fig. 9
figure 9

The yielded elements and stress histories after producing both the lower and upper trona seams: woDESS (a), RN_30 (b), RN_60 (c), RN_90 (d), stress history in the upper (T3) (e) and lower (T4) stiff layers when the upper stiff layer is destressed

In the second stage of the numerical study, alternative borehole patterns were investigated. Both BH_STR and BH_STG patterns in the upper stiff layer (T3) increased the stress in between the weakened rock columns, similar to RN_30 and RN_60 (Fig. 10). Compared to the model with no destressing (6.6 MPa), 22% stress increase in BH_STG, and 46% stress increase in BH_STR were observed where the maximum principal stresses were 8.0 MPa and 9.6 MPa respectively. Unlike the ribbon pattern, no stress relaxation was observed in the lower stiff layer (T4) (Fig. 10). For BH_STR case, 19% stress increase in lower stiff layer (T4) and similar stress values in BH_STG case were observed relative to the model simulated without destressing (woDESS case). Simulations point out that both cases are successful in destressing that drives the failure of the stiff stratum due to the increasing stress accumulation.

Fig. 10
figure 10

The yielded elements when BH_STG applied and after the first production stage (a), producing the whole deposit (c); when BH_STR applied and after the first production stage (b), producing the whole deposit (d), stress history in the upper (T3) (e) and lower (T4) (f) stiff layers with different borehole patterns

To conclude, both BH_STG and BH_STR patterns were useful for destressing. In between the destressed rock columns, a stress increase by 22% and 46% implies that the yielded zone around the boreholes induce the propagation of fractures through the rest of the rock mass. Although BH_STR shows a higher stress accumulation, BH_STG pattern also works. By considering the operational efficiency the BH_STG pattern can be concluded to be more suitable for destressing.

Assessment of the simulations with regard to the economic and operational considerations shows that the ribbon pattern with 60 m off-set (RN_60) and the staggered borehole pattern (BH_STG) are favorable for destressing the upper stiff layer (T3). Further, these patterns were compared to determine the best method for the elimination of sudden collapse risk. The stress history shows that BH_STG is superior compared to RN_60 regarding the stress increase in upper stiff layer (T3). The maximum principal stresses are 7.2 MPa for RN_60 and 8 MPa for BH_STG with respect to 6.6 MPa in the model simulated without destressing (woDESS case). Although both patterns are suitable for inducing fracture propagation due to stress accumulation (10% in RN_60, 22% in BH_STG), BH_STG pattern is a better option in the stiff and thick T3 layer considering time & cost constraint, and the operational advantages.

Besides the effects of pattern size and shape, the effect of destressed layer was investigated for multiple stiff layers. The model was designed to involve two distinct strata with stiff mechanical properties to determine the best layer for performing the destressing operation. Models were arranged to simulate the best practice for both of the patterns. Similar to the destressing in upper stiff layer (T3), BH_STG and RN_60 patterns were implemented to the lower stiff layer (T4). The pattern dimensions were adjusted regarding the thickness of the T4 layer while the horizontal plane dimensions were identical. A similar DFN model and joint mechanical properties were used.

Investigating the stress history in the upper stiff layer (T3) (Fig. 11) for RN_60 pattern, σ1 values increase by 10% due to the stress transfer from the weakened zone. However, if the lower stiff layer (T4) was destressed, it was observed that σ1 values in the upper stiff layer (T3) decreases by 8%. In the lower stiff layer (T4), stress relaxation is around 21% due to the implementation of DFN into the upper stiff layer (T3). However, destressing the lower stiff layer (T4) induces only a 3% stress increase. This trend is similar to the major principal stress values in T3 when destressing implemented. However, in the former case, the stress increase was 10% with respect to the 3%. The results show that both cases are acceptable since the strata fails in a controlled manner in either of the cases.

Fig. 11
figure 11

The major principal stresses in T3 and T4 layers for destressing with T3_RN_60 and T4_RN_60 patterns

For BH_STG pattern, destressing the upper stiff layer (T3) induces a 22% increase in the stress accumulation of T3 layer; however, it is a 7% relaxation when the lower stiff layer (T4) is destressed (Fig. 12). These outcomes similar with the RN patterns as implementation of destressing into the upper layer results in higher stress levels in the respective layer and lower stress levels in the lower layer. However, σ1 as a result of destressing the lower stiff layer (T4) shows 2% relaxation when BH_STG is conducted in T4 and no significant σ1 change is observed in T3.

Fig. 12
figure 12

The major principal stress change in the upper (T3) and lower (T4) stiff layers with respect to the destressing pattern and location

Figure 13 shows the relaxation in stiff layers achieved by BH pattern. The difference between the major principal stress levels in RN and BH patterns can be explained by the geometry and size of the weakened volume. The stress relaxation in RN patterns is sharper than BH patterns but both patterns are successful in preventing the sudden collapse risk. Regarding the simulations, for deep deposits it is more convenient to conduct destressing operation within the stiffer and thicker layers such as it is done in the upper stiff layer (T3) in this case. By increasing the stress accumulation in the unweakened zone, the failure of the strata is ensured.

Fig. 13
figure 13

The major principal stresses in T3 and T4 layers for destressing with T3_BH_STG and T4_BH_STG

Table 4 provides a summary of the numerical outcomes, which covers;

  1. 1.

    determination of the best ribbon and borehole patterns, when the upper stiff layer is destressed.

  2. 2.

    stress relaxation or accumulation on the lower stiff layer, when the upper stiff layer is destressed.

  3. 3.

    comparison of the best ribbon and borehole patterns.

  4. 4.

    manipulation of the stress field on the undisturbed stiff layer when the other stiff layer is destressed.

Table 4 Summary of the research outcomes

Conclusion

This study explores alternative destress blasting patterns to mitigate the sudden collapse risk in solution mines with stiff roof strata. Parametric simulations of ribbon and borehole destressing patterns were implemented in a DEM code with calibrated geomechanical inputs derived from field and laboratory data. Simulations reveal that the sequential removal of trona layers leads to stress accumulation in stiff roof layers, particularly T3 and T4, while significant vertical displacements in lower layers (T5 and T6) imply bending and failure. Asynchronous deformations in the cavern roof strata poses a sudden collapse risk. The best stress relaxation was observed in implementing 60 m evenly spaced ribbons within the upper stiff layer, which reduced stress in the lower stiff layer by 21%. On the other hand, the staggered borehole pattern provides the necessary stress redistribution to induce fracture growth by increasing stress concentration both on the upper and lower stiff layers. As both patterns provide effective stress management for directly or indirectly inducing stress relaxation, operational and economic concerns should be regarded to decide on the feasible method. Setting destress panels with regular off-set distances requires creating an intensely fractured rock mass on large scale. However, borehole staggered pattern may induce a similar fracture network by increasing stress concentration in-between fragmented rock columns. Another research highlight is the significance of destressing location. Destressing in the upper stiff layer (T3) reduced stress concentration by 24% and facilitated stress relaxation across multiple layers, highlighting its strategic importance.

Potential future research directions are exploration of multi-cavern interactions, incorporating advanced numerical modeling techniques with hydromechanical or thermomechanical coupling, and validating findings with field-scale monitoring to further enhance the reliability and applicability of the proposed destressing methods. Expanding destress research with regard to these aspects will enable more comprehensive risk mitigation and improve the long-term sustainability of solution mining operations.