*Article* **A Field Study on the Law of Spatiotemporal Development of Rock Movement of Under-Sea Mining, Shandong, China**

**Jia Liu 1,2,3, Fengshan Ma 1,2,\*, Jie Guo 1,2, Guang Li 1,2, Yewei Song 1,2,3 and Yang Wan 1,2,3**


**Abstract:** The phenomenon of rock movement in mining areas has always been a difficult problem in mining engineering, especially under complicated geological conditions. Although the backfilling method mitigates the destruction of the surrounding rock, deformation can still exist in the mining area. In order to ensure the safety of under-sea mining, it is necessary to study the rock movement laws and the mechanisms. This paper focuses on a settlement analysis of the monitoring data of the No. 55 prospecting profile. By analyzing the shape of the settlement curves, the spatial distribution characteristics of settlements of different mining sublevels are summarized. Additionally, the fractal characteristics of the settlement rate under different space–time conditions are studied. We also discuss the relationship between the fractal phenomenon and the self-organized criticality (SOC) theory. The findings are of great theoretical value for the further study of mining settlements to better understand the physical mechanisms of internal movement and rock mass failures through the fractal law of the settlement. Furthermore, elucidating the rock movement law is an urgent task for the safety of seabed mining.

**Keywords:** rock movement; high-dip metal deposits; fractal character; the self-organized criticality (SOC) theory

#### **1. Introduction**

Rock movement of different degrees often occurs in mining areas due to the disturbance of mining activities [1–4]. The movement, deformation, and destruction of rock masses caused by underground mining activities may further lead to ground subsidence which has caused growing social stability and ecological problems worldwide [5–7]. Therefore, many scholars have been attaching great importance to the study of rock movement phenomena; they have proposed relevant theories and numerical models [8–10]. They have conducted much research on rock movement law and relevant profound theories from different perspectives [8,11–13]: geometric theory, mechanical theory, and mathematical–statistical theory, for instance, the geometric method [14], the circular integral grid method [15], key strata theory [16,17], random medium theory [18], numerical methods [19,20], and fractal theory [21,22]. In addition, a large number of numerical models have emerged, including the Knothe function model [23], the visco-elastic model [24], the anisotropic material model [25], FDM predictive methodology of subsidence [26], the improved Knothe function [27–29], the "four-zone" model [30], and the improved key stratum theory [31]. Regarding fractal theory, the fractal concept was created in 1967 [32]. For the study object, if there is a power-law (or scaling-law) relationship between some of its variables within a certain range, it is considered to have a fractal phenomenon [33,34]. The power-law and fractal theory are closely linked and complement each other. The fractal theory has been widely applied in social sciences [35], agriculture [36], geography [32],

**Citation:** Liu, J.; Ma, F.; Guo, J.; Li, G.; Song, Y.; Wan, Y. A Field Study on the Law of Spatiotemporal Development of Rock Movement of Under-Sea Mining, Shandong, China. *Sustainability* **2022**, *14*, 5864. https:// doi.org/10.3390/su14105864

Academic Editors: Longjun Dong, Yanlin Zhao and Wenxue Chen

Received: 17 March 2022 Accepted: 11 May 2022 Published: 12 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

geology [37,38], and other fields in recent years. Although the application of fractal theory in geological science started relatively late, it shows strong vitality and broad application prospects due to its powerful expression and concise form. Most geological phenomena are fractal in the statistical sense. In addition, the research and exploration of rock mechanics have also benefited from the development of fractal theory [22,35,36], for instance, in the study of the fractal characteristics of land subsidence [16] and the acoustic emissions of rock mass [37–39]. Fractal theory is widely used to guide production practices.

In summary, rock movement theory continues to develop and mature. However, most mining-induced rock movement studies focus on coal mines; the problem of rock movement caused by the mining of high-dip metal deposits on the seabed is little studied, especially in the presence of faults [40–44]. As one of the popular mining methods, backfill mining is increasingly applied in China and other mining countries. The filled materials consolidate and compress together with the surrounding rock, and the stress of the surrounding rock gradually shifts to the filling body over time. Finally, the mine reaches a state of equilibrium. However, the movement of the surrounding rock is inevitable, which is affected by the mining operation, the low-stiffness filling body, the filling method, the existence of unfilled void spaces, and repeated mining activities [10,45,46]. Severe rock movement or even surrounding rock damage can pose a great threat to the safety of mining personnel and mining production activities. Nonlinear scientific theories are increasingly widely used in geological hazard evaluation and prediction [47–51]. Some scholars have used the self-organized criticality (SOC) theory to analyze the progressive failure process of rocks under stress. We know that many natural phenomena exhibit fractal features, and geological phenomena are no exception. As for the physical mechanisms underlying the fractal phenomenon, the self-organized criticality (SOC) theory offers a better explanation of the characteristics of nonlinear complex systems [47,48,51]. Simply put, SOC means that a system will spontaneously evolve to a critical state. The stratum is in a dynamic equilibrium state before excavation. The excavation process breaks the equilibrium and leads to changes in the energy, so the stratum instinctively adjusts itself and organizes itself to reach a new equilibrium through the adjustment of stress and deformation.

In this study, we used the following methods to explore the law of rock movement produced by mining high-dip seabed deposits. We took the No. 55 prospecting line of the Xinli mining area as the study object, summarized the deformation characteristics of various sublevels (−200 m sublevel, −240 m sublevel, −320 m sublevel, −400 m sublevel, −480 m sublevel, and −600 m sublevel), analyzed the vertical settlement rate using fractal theory, and studied the possible deformation mechanism using the self-organized criticality (SOC) theory. Hereafter, further studies on the mechanisms of rock movement phenomena caused by mining high-dip deposits will be conducted, including physical simulation and numerical simulation research.

#### **2. Study Area**

#### *2.1. Geographical Environment and Geology*

The Xinli mining area of the Sanshandao gold deposit is the first hard rock mine that belongs to the under-sea deposit in Shandong province, China. It is located in Laizhou bay, Laizhou city, bordering the Bohai Sea in the northwest and connected to the land in the southeast (Figure 1). Many scholars have conducted much research on it [52–55]. The study area belongs to the continental semi-humid climate of the temperate monsoon region. In addition, it has both oceanic and continental climate characteristics. The annual average temperature is 12.5 °C, and the precipitation is 612.1 mm [56].

**Figure 1.** The study area (indicated in the red circle).

The study area is located to the east of the Tan–Lu fault, the north-limb of the Qixia anticlinorium, and the northeastern part of the Sanshandao–Cangshang fault zone. As for the Xinli mining area, the gold-bearing alteration zone is 1300 km long in the strike, and its width varies from 70 to 185 m. It is produced in a relatively simple large vein shape, with a stable occurrence and good mineralization continuity. The altered rocks are zoned along the strike and inclination. From top to bottom, the altered rocks are pyrite sericite granite, pyrite sericite granitic cataclysm, fault gouge (existing below the main fracture), pyrite sericite cataclysm, pyrite sericite granitic cataclysm, and pyrite sericite granite. The pyrite-sericite cataclysm zone in the range of 0–35 m below the main fracture has the strongest alteration and gold mineralization, and it is also the location of the main ore body. The quaternary sediments are distributed in the shallow part of the surface, exposed in the southeast, with a thickness of 16–46 m, with sandy clay, fine sand, medium-coarse sand, and coarse gravel. The sedimentary sequence is diverse, and the regularity is poor [57]. The direction of the maximum principal stress is roughly the horizontal direction (NWW), and its value is greater than the self-weight stress of the rock mass [7]. There are three main faults (F1, F2, and F3) and some typical secondary joints in the study area (Figure 2). The gold ore body is an altered rock type gold deposit controlled by the F1 fault (a reverse fault), whose strike is 40◦ and the dip angle is about 38◦. The F2 fault (a normal fault) lies west of the F1 fault, with a strike of 15◦ and a dip angle of 80◦. The F3 fault (a strike-slip fault) stretches northwest to the sea and southeast to the continent, with a strike of 300◦ and a dip angle of 90◦.

**Figure 2.** Three main faults and the Xinli mining area.

#### *2.2. Mining Situation*

The Xinli mining area of the Sanshandao gold deposit began its basic construction in 2005, with a designed production capacity of 1500 t/d. The Xinli mining area adopted the mechanical upward horizontal slice stopping-filling method with pointed pillars. In order to improve production efficiency, simultaneous upward mining in several sublevels has been designed. The mining activities are conducted from bottom to top by slicing and filling at a slice height of 2.5 m. After the excavation, the mined-out area is immediately backfilled with hardened solids mixed with mineral tailings and cement slurry (binder). The stope length is 100 m, and its width is the thickness of the ore body. Three vertical shafts and one return air shaft have been established in the mining area, and several sublevels (−200 m, −240 m, −320 m, −400 m, −480 m, and −600 m.) have been developed (Figure 3). In addition, a transportation roadway was developed at −440 m; a water silo, pump room, and power distribution room at −600 m; and a fine ore recovery roadway at −667 m.

**Figure 3.** The simplified geological condition of the No. 55 prospecting profile (the green lines show the plot range of the following settlement curves).

#### **3. Methodology**

#### *3.1. Displacement Monitoring of the Roadway*

In August 2015, a settlement monitoring system was established in the Xinli mining area. In the whole mine, only the No. 55 prospecting profile line passes through the ore body vertically from the hanging wall, that is, perpendicular to the mining direction, forming the most favorable monitoring space conditions (Figure 2). Compared with other profiles, the monitoring network of the No. 55 prospecting profile is relatively complete, covering several sublevels of different depths. The monitoring network is mainly deployed in several sublevels of the No. 55 prospecting profile and other roadways, and the number of monitoring points exceeds 250 (Figure 4); the monitoring intervals are usually one or two months. All of the monitoring uses the fourth-level leveling method. A leveling survey is a method used to determine the height difference between different points on the ground by using instruments (such as level and leveling staff) that provide level realization. It is a vital technique of classical geodetic surveying [58–60]. This method played a positive role in China's height measurement campaign of Mount Qomolangma from 2019 to 2020 [61]. The total length of the leveling route is 7100 m. In the early stage, the monitoring interval was one month. Later, due to the increased monitoring workload and relatively stable settlement data in the study region, the monitoring interval was adjusted to two months after 2017. In order to ensure the efficiency of the monitoring work and to save costs, the distance between the monitoring points was nearly 50 m, and it became smaller in the area close to the ore body, where the deformation is intense. The density of the monitoring points depends on the deformation degree of surrounding rock; areas of intense deformation often require more detailed monitoring data. It should be noted that the deformation of the roadway caused by the mining disturbance is composed of the convergence displacement

of the roadway itself and the displacement of the rock mass. When considering the time effect, it can be assumed that the convergence displacement of the roadway itself is close to zero. Therefore, the monitoring data of the roadway were taken as the results of the rock movement.

**Figure 4.** The monitoring network of the No. 55 prospecting profile.

#### *3.2. Fractal Theory*

Regarding fractal theory, the capacity dimension is also called the box-counting dimension, and its definition is based on coverage [62]. The functional model is expressed as N (r) = C × <sup>r</sup><sup>−</sup>D. In the equation, r represents characteristic scale; C represents the proportional factor; D represents the fractal dimension; and N (r) represents the number or sum of the scale greater than or equal to r. According to the different meanings of N (r), it can be divided into two function models of the capacity dimension: the content–number method and the content-sum method [63]. In order to obtain the fractal dimension value D, the value of r is constantly changed, and a series of corresponding N (r) values are obtained. The two series of data obtained are projected onto a logarithmic coordinate graph. If there is a straight line in a log–log plot, it indicates that the data conform to the fractal distribution. In order to analyze the fractal characteristics of settlement rate evolution in the Xinli mining area, we adopted long-term monitoring data of the −200 m, −320 m, −400 m, and −480 m sublevels from May 2016 to November 2021.

#### **4. Results**

#### *4.1. Settlement Monitoring Curve*

Although it is well known that the backfilling mining method is an effective method used to control the severe deformation of surrounding rock [64–67], the rock movement phenomenon still exists in the mining region. Furthermore, the strength of the surrounding rock has decreased due to mining disturbance. We reorganized the monitoring results of several sublevels of the No. 55 prospecting profile from May 2016 to November 2021. The monitoring curves are displayed in Figures 5–7. It should be noted that different colors imply different monitoring dates, and we only present some of the monitoring dates in order to keep the figures more concise and readable. Due to the arrangement of the administrative staff, the monitoring of the −240 m sublevel was interrupted in April 2019, and the monitoring of the −600 m sublevel was started in December 2017. In addition, the settlement value increased by 1000 times in the following figures, and the horizontal distance is the real value. From the final form of the settlement monitoring curves, the backfilling method cannot completely constrict rock movement and deformation of the surrounding rock. According to the years of the monitoring results, we found that

the vertical settlement curves can be divided into two categories based on the shape of the curves: (1) a "pan"-type curve. The upper sublevels (−200 m and −240 m) of the No. 55 prospecting profile tend to show "pan" type curves (Figure 5). Most monitoring points have similar settlement values, and the curves are smooth. (2) a "funnel"-type curve. The lower sublevels (−320 m, −400 m, −480 m, −600 m) of the No. 55 prospecting profile tends to show "funnel"-type curves (Figures 6 and 7). The settlement values of most of the monitoring points have obvious differences, and the curves are somewhat rough. Although there are also some differences between the two curves, there are some similarities, regardless of the sublevels of the No. 55 prospecting profile; the monitoring points on the left side of the monitoring line have a slight upward trend. After mining, the surrounding rock moves and deforms as a result of disturbance, and the closer the monitoring point is to the mining area, the larger the displacement is. Moreover, the maximum settlements of the six monitoring sublevels from top to bottom are 296 mm, 146 mm, 308.5 mm, 332.5 mm, 268.5 mm, and 82.5 mm, respectively. In the −200 m, −320 m, −400 m, and −480 m sublevels, the horizontal distances between the monitoring points reaching a 150 mm settlement are 341.6 m, 251.3 m, 152.3 m, and 97.9 m, respectively.

**Figure 5.** The "pan"-type curve (**a**): −200 m sublevel (**b**): −240 m sublevel.

**Figure 6.** The "funnel" type curve (**a**): −320 m sublevel (**b**): −400 m sublevel.

**Figure 7.** The "funnel" type curve (**a**): −480 m sublevel (**b**): −600 m sublevel.

#### *4.2. Fractal Character of the Settlement Rate*

Studies have shown that the dynamic settlement curve of the ground surface points caused by mining has statistical fractal characteristics [22,68].

We used two function models of the capacity dimension to analyze the monitoring data (Tables 1 and 2). The value of r (settlement rate) was constantly changed, and a series of corresponding N (r) values were obtained. It should be noted that the absolute values of the settlement are used here for convenience.

**Table 1.** Content-number method.


**Table 2.** Content-sum method.


The two series of data obtained were projected onto a logarithmic coordinate graph, respectively (Figure 8). For the content-number method, there is a straight line that was fitted by the least square method in a log–log plot, with an R<sup>2</sup> value of 0.95. For the contentsum method, there is a straight line that was fitted by the least square method in a log–log plot, with an R<sup>2</sup> value of 0.92. Whichever method we used, they both indicated that the data conformed to the fractal distribution.

The settlement rate data of any provided monitoring sublevel from May 2016 to November 2021 also showed a fractal character. When taking the −400 m sublevel as an example, the settlement rate data also have an obvious fractal character (Table 3 and Figure 9). There is a straight line that was fitted by the least square method in a log–log plot, with an R2 value of 0.97. In addition, the settlement rate data of the −400 m sublevel from June 2016 to September 2019 also showed a fractal character (Table 4 and Figure 10). The settlement events over different space–time scales both show fractal characteristics. The settlement events reflect the scale invariance of the settlement system.

**Figure 8.** The settlement rate distribution of the whole No. 55 prospecting profile from May 2016 to November 2021. (**a**) Content-number method; (**b**) content-sum method.

**Table 3.** The settlement rate distribution of the −400 m sublevel from May 2016 to November 2021 (content-number method).


**Figure 9.** The log–log plot of the settlement rate distribution of the −400 m sublevel from May 2016 to November 2021.

**Table 4.** The settlement rate distribution of the −400 m sublevel from June 2016 to September 2019 (content-number method).

**Figure 10.** The log–log plot of the settlement rate distribution of the −400 m sublevel from June 2016 to September 2019.

#### **5. Discussion**

The subsidence monitoring curves at different depths have distinct characteristics. As the depth increases, the difference in settlement values between the monitoring points becomes larger. As mentioned earlier, the horizontal distance between the monitoring points reaching 150 mm subsidence is decreasing. There are many factors that could contribute to this result: the mining depth, the thickness of the ore body, the mining history, the mining method, backfilling effect, and tectonic geological conditions. Furthermore, the spatial accumulative effect may be another factor. Tectonic stress, hydrogeological conditions, and temperature vary at different mining depths; the thickness of the ore body determines the excavation scale; the mining history determines the original disturbance and deformation of the surrounding rock; the mining method and backfilling method also affect the settlement of the surrounding rock; and ore-controlling fractures have a significant influence on the deformation of mining areas [41]. All of these factors could lead to the result that the upper part may suffer more severe disturbances. Dynamic disturbances often cause local instability of the surrounding rock and goaf, which further lead to different degrees of settlement. When the lower surrounding rock is disturbed by mining, deformation occurs, and the filling is not sufficient to form goaf. The upper part tends to deform by the influence of the excavation of the lower part. This is a natural occurrence that happens under the influence of gravity. Moreover, the upper part can deform when mining activities are implemented in the upper rock mass. This means that the uppermost sublevel is the most affected by the spatial accumulative effect.

Based on the results of settlement monitoring curves, we also found apparent asymmetric characteristics, similar to other scholars who have conducted much research on subsidence [69–72]. The steep geometry of the orebody and the existence of the fault both have potential effects on rock movement phenomena. The steep geometry of an orebody and the existence of a fault can affect rock movement by forming different stress conditions. In this paper, the monitoring points on the right side of the monitoring line have a downward trend; the monitoring points on the left side of the monitoring line have a slight upward trend. We hold that the flexural cantilever beam model [73,74] is suitable to explain this phenomenon. Figure 11 is the schematic diagram. At first, the inclined orebody and the surrounding rock are under prominent horizontal tectonic stress [75]. Due to mining disturbance and tectonic stress, the deep rock mass produces a trend of upward bending, and the shallow rock is pulled; then, the fault activation of the shallow surrounding rock under tensile stress leads to the formation of a sedimentary basin. The hanging wall and footwall fault blocks can be regarded as two bending cantilevers interacting with each other. The equilibrium force generated by extension causes cantilever uplift of the footwall fault block and settlement of the footwall fault block.

**Figure 11.** Schematic diagram of settlement development.

The rock mass is heterogeneous and consists of rock blocks and structural planes. Because of the difference in the physical and mechanical properties of rock block, as well as the complexity of the structural plane, mining-induced rock mass failures are strongly nonlinear, and the dynamic settlement curve of the surrounding rock is fluctuant. Taking the "J12" monitoring point in the −400 m sublevel as an example, the settlement rates show obvious fluctuation characteristics. (Figure 12). Mining-induced rock mass failures show a strong nonlinearity, which makes the dynamic settlement curves of the No. 55 prospecting profile show a non-smooth shape (Figures 5–7). The evolution process of rock movement is essentially a self-organizing evolution process, from disorder to order and from low levels to high levels under the action of certain external conditions and internal nonlinear mechanisms [49]. Underground mining activities may cause movement of the surrounding rock due to the advance of the excavation fronts and the progressive closure or collapse of the mineral extraction galleries. During this excavation process, the stress balance of the rock mass was disturbed, and stress redistribution occurred. The rock movement system of the mining area is an open system because the exchange of material and energy is always occurring between the rock movement system and the outside. The surrounding rock stress field, rock mass structure, and other factors all change with respect to two scales: time and space; and the nonlinear interaction between them makes this system very complicated. Thus, the rock movement system remains far from being an equilibrium, although it may have stable statistical properties [75]. The fractal character of subsidence is related to the SOC of the rock movement system. Thus, the internal movement and failure mechanisms of the rock mass are still unexplained problems. This paper reveals the fractal characteristics of rock movement, similar to other natural phenomena. According to the fractal characteristics and the SOC theory, we may find rock movement anomalies in time. In addition, it is necessary to use fractal theory to conduct in-depth studies on rock movement mechanisms.

**Figure 12.** The settlement rates of the "J12" monitoring point in the −400 m sublevel.

#### **6. Conclusions**

We used several years' worth of monitoring data of the No. 55 prospecting profile of the Xinli mining area to study the deformation characteristics. The results show that the dynamic settlement curves of the monitoring sublevels have interesting differences. The vertical settlement curves can be divided into two categories based on the shape of the curves: (1) a "pan"-type curve (−200 m and −240 m sublevels) and (2) a "funnel"-type curve (−320 m, −400 m, −480 m, and −600 m sublevels). This phenomenon is not only related to the basic geological and mining conditions but also to the cumulative effect. The upper part may deform as a result of the mining activities in this part as well as due to the deformation of lower sublevels. This is a natural occurrence that happens under the influence of gravity.

Mining-induced rock mass failures are strongly nonlinear due to the randomness of rock mass structure, which is reflected in the fluctuating settlement curves. In addition, the settlement rates of the monitoring points of the No. 55 prospecting profile have good statistical fractal characteristics. The settlement events at different spatial and temporal scales all have fractal characteristics, which reflects the scale invariance of the settlement system. Based on the above information, the self-organized criticality (SOC) theory offers a better explanation of the fractal characteristics of nonlinear complex systems. In this paper, the physical mechanisms of rock movement have been preliminarily investigated through the collection and analysis of the monitoring data, and the findings have theoretical value for the further study of mining settlements.

**Author Contributions:** Conceptualization, J.L. and F.M.; methodology, J.L. and J.G.; formal analysis, J.L. and F.M.; investigation, J.G. and G.L.; resources, Y.W.; data curation, Y.S.; writing—original draft preparation, J.L.; writing—review and editing, F.M.; supervision, J.G.; project administration, G.L. and Y.W.; funding acquisition, F.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by The National Natural Science Foundation of China (Grants number 41831293, 42072305).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to express their sincere gratitude to the Sanshandao Gold Mine for their data support. In addition, the authors are grateful to the assigned editor and the anonymous reviewers for their enthusiastic help and valuable comments, which have greatly improved this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

