Owing to the development of coal mining technology and mechanical manufacturing in China, high-strength mining has become the mainstream modern-mining method. Extensive mining damage is accompanied by severe overlying fractures and ground-fissure evolution (Li et al. 2022). Because the strength of the loose layer (LSL) is significantly lower than that of the overlying rock (OLR), the mining damage on the LSL and OLR differs significantly. It primarily manifests as a continuous flow until it stabilizes with the movement of the OLR, which renders it difficult to predict the evolution and development of ground fissures.

To study the law of surface movement and effectively control the phenomenon of surface failure under the action of high-strength mining damage, many scholars regard the occurrence characteristics of sand particles with strong fluidity as random media and analyze their mutual response motion with OLR breakage. According to the evolution law of surface damage under the influence of high-strength mining in super large longwall faces, it is found that surface deformation, ground fissures, and soil physical properties in aeolian sand areas have certain self-healing characteristics (Hu et al. 2016; Xu et al. 2017; Wang et al. 2015). Several researchers have analyzed the sensitivity and deformation characteristics of mining influence on the surface movement and deformation of the windy and sandy region based on on-site measurements (Hu et al. 2016; Wang et al. 2024; Xu et al. 2016; Yang et al. 2019) and reveal the evolution law of dynamic movement parameters in subsidence basin and ground fissure (Liu 2023; Yang et al. 2023; He et al.2021). Besides, due to the compression characteristics and fluidity between particles in the loose layer, the ratio of soil layer to bedrock has a significant impact on surface movement (Liu et al. 2022; Zhang et al. 2016). As the ratio continues to increase, the surface deformation shows a clear characteristic that increases first and then decreasing. However, when it reaches a certain limit, the surface deformation tends to stabilize (Liu et al. 2022; Hu et al. 2015).

In response to the particularity of asymmetric surface subsidence in extremely high mining faces, scholars have improved the probability integration method for subsidence prediction based on the control effect of thick-hard rock layers on surface subsidence. They believe that increasing the LSL thickness intensifies the asymmetry of surface subsidence (Gao et al. 2018; Yu et al. 2018). Therefore, the variation characteristics of the LSL during the mining are key factors affecting the surface movement law. To systematically analyze the variation characteristics of the control factors of the OLR and the passive factors of the LSL during damage conduction, researchers have investigated the stress evolution law under the bearing capacity of loose arch structures. The dynamic relationship between the primary key stratum (PKS) movement and loose arch-bearing capacity was revealed (Wang et al. 2019; Sun et al. 2019 and 2020). The effects of parameters such as the lateral pressure coefficient, internal friction angle, and cohesion on the shape characteristics of loose arch structures were investigated via numerical simulations (Wang et al. 2019). To address the differences in the movement of the OLR-LSL double medium under the effect of mining, the effects of the OLR movement state and LSL characteristics on surface movement were revealed base on a subsidence model under the composite action of the OLR-LSL and a prediction formula (Hu et al. 2022; Liu et al. 2019, 2023; Yang et al. 2023; Gu et al. 2013; Luo et al.2023). To address mining damage conduction from the longwall face to the surface, scholars constructed a “space-sky-surface” integrated multiple asynchronous monitoring method (Banerjee and Raval 2022; Simmons and Wempen 2021; Zhang et al. 2021; Cai et al. 2023) and proposed a “3-3-1” relationship among the main roof breakage distance, the PKS breakage distance, and the distribution distance of main ground fissures, and explained the influence of advancing speed on OLR-LSL deformation (Li et al. 2021; Wang et al. 2024).

Owing to the concealment of geological conditions and the randomness of mining-damage transmission in the LSL and OLR, on-site monitoring cannot capture the development of mining fractures and the characteristics of rock fractures and collapse in a timely manner owing to the limitations of its monitoring range and period. Therefore, numerical simulation methods are widely used to investigate OLR fracture mechanisms, soil movement, and deformation characteristics owing to their convenience. Currently, numerical simulations are primarily based on continuum methods (including the boundary-element method (BEM), finite element method (FEM), and Lagrangian finite-difference method (FDM)), as well as discontinue methods (including the discrete-element method (DEM) and discontinuous deformation analysis (DDA)) and the coupling methods of FEM/DEM and FDM/DEM) (Lisjak and Grasselli 2014; Gao et al. 2014; Huang et al. 2023; Yan et al. 2023; Breugnot et al. 2016; Yin et al. 2020). the coupled numerical simulation methods of FEM/DEM and FDM/DEM combines the advantages of FEM and DEM, but it requires finite element mesh division for each block, and the computation consumption is extremely high, especially for 3D simulation. It is difficult to apply to engineering scale simulation, especially for engineering problems with strong disturbances such as coal mining. Thus, studies pertaining to the microconduction mechanism of mining damage in the OLR-LSL are primarily performed using the DEM via the Universal Distinct Element Code (UDEC) and Particle Flow Code (PFC) software (Gao et al 2014 and 2022; Rorato et al. 2021; Yan et al. 2023; Zhang et al. 2023; Qin et al. 2023). The OLR and LSL are generally modeled using block and particle elements, respectively. However, researchers disregarded the distinct characteristics of two discontinuous media and regarded them as the same discrete-element medium for the simulation analysis of OLR structural instability, fracture evolution, and surface subsidence (Cheng et al. 2020; Wang et al. 2018; Wei et al. 2022; Xu et al. 2023; Lou et al. 2021). This can result in inaccurate descriptions of the mining damage and development characteristics (Gao et al. 2014).

Thus, based on the respective characteristics of OLR and the LSL, a combining simulation method based on block and particle discrete elements was devised in this paper to analyze the collaborative movement characteristics of OLR and the LSL. The conduction characteristics of mining damage under a double medium are simulated and analyzed to determine the effect of the OLR breakage morphology on surface subsidence and to reveal the mining-damage conduction mechanism.

The overburden in coal mining includes two primary types of media: the OLR and LSL. The overall relative strength of the OLR is relatively high, and because of its internal joint-layered structure, it propagates in the form of a combined collaborative body under the control of a high-stiffness rock layer (the PKS). The LSL contains scattered particles composed of three-phase materials, with almost zero tensile stress and cohesion between the particles. It exhibits strong fluidity, and uneven displacement occurs between the particles after the PKS becomes unstable (Liu et al. 2022; Gu et al. 2013; Gao et al. 2014). Previous studies showed that the DEM is an effective method for simulating the physical deformation of rocks and granular materials and can effectively simulate the large deformation characteristics of mining activities (Chen and Hu 2018; Sun et al. 2019 and 2020; Wang 2018; Khanal et al. 2017). In the UDEC simulation method, although rigid blocks can effectively simulate the fracture of OLR and the development characteristics of mining-induced fractures, they can only slide along a predetermined fracture surface (Gao et al. 2014). Thus, it is not suitable for simulating the movement and deformation characteristics of LSL particles. In the PFC simulation method, although the particle model can reproduce the development characteristics of the OLR fractures well, the simulation effect of the rock fracture structure is unsatisfactory. In addition, because of the limitation of the model size, if the model is fully modeled under a certain particle size, then the computational efficiency will be degraded (Zhao and Wang 2011; Feng and Owen 2014). Therefore, considering the interaction between the OLR and LSL, the OLR determines the bottom boundary of the LSL, and the LSL determines the top load of the OLR. A simulation analysis was performed on the combining movement results of the OLR-LSL. The modeling and data-transfer methods between the block and particle discrete-element simulation methods are illustrated in Fig. 1.

The modeling and data transfer methods between UDEC and PFC^2D simulation method

In terms of the computational mechanics of rock mass, the DEM can accurately reflect geometric features (Cheng et al. 2020; Sun et al. 2019; Gao et al. 2014). Furthermore, it is suitable for addressing rock mass failure and large deformation problems caused by nonlinear deformation and concentrated failure at the contact between blocks. In this study, the block DEM software, UDEC (Itasca Consulting Group Inc. 2013) was adopted to investigate the characteristics of OLR movement caused by longwall coal mining activities, whereas the particle DEM software, PFC, was used to investigate the characteristics of LSL surface subsidence.

The fracture and constitutive behavior between blocks for the block DEM are shown in Fig. 2. When the normal tensile stress is greater than the tensile strength of the joint t, the joint will undergo tensile failure and the normal stress σ_n = 0. The shear stress τ_s is determined by the cohesion c and friction angle φ between the contact blocks (Eq. (1)). The deformation of the joints is determined by contact stiffnesses K_n and K_s. When the shear stress applied to a joint exceeds the maximum shear strength, the joint undergoes shear slip failure (Itasca Consulting Group Inc. 2013; Lisjak and Grasselli 2014).

where \(\Delta \mathop \mu \nolimits_{\text{s}}^{\text{e}}\) is the elastic region of the shear-displacement increment, K_s is the shear stiffness, c represents cohesion, and φ is the friction angle.

UDEC modeling of fracture development in strata. a Normal and shear stiffnesses between blocks; b Constitutive behavior in shear and tension

Because the LSL in the study area is primarily composed of aeolian sand, the bonding strength between particles is negligible, and the effect of particle shape is primarily considered at the contact between the particles (Wensrich and Katterfeld 2012; Yin et al. 2020; Rorato et al. 2021). Thus, the rolling resistance contact model was used in the PFC simulation, which included a conventional linear elastic–frictional contact model for particle relative displacement at the contact and an additional set of elastic spring nontensional joints and sliders for the rolling motion (Fig. 3). The contact stiffnesses are expressed as follows (Jiang et al. 2005; Itasca Consulting Group Inc. 2012b):

where E_mod and K_ratio are the material parameters to be calibrated, D the diameter of the smallest contacting sphere, L the distance between the grain centers, K_n the normal stiffness, and K_s the shear stiffness.

where K_r is the rolling stiffness, R the effective radius; and R₁ and R₂ the radii of the two particles in contact.

where μ_r is the rolling friction coefficient and F_n is the normal contact force.

PFC modeling of rolling resistance contact model. a Normal, shear and rolling stiffnesses between particles; b Elastic-perfectly plastic model accounting for rolling resistance at contact

The mechanical interaction and combining between the OLR simulation in the UDEC and the LSL simulation in the PFC included the load of the LSL on the OLR and the effect of the OLR movement on the LSL. Thus, the bottom stress from the LSL simulation (using the PFC) and the top breakage morphology from the OLR simulation (using the UDEC) must be extracted, as shown in Fig. 1. First, the model was constructed based on actual geological conditions and included two components: the PFC construction of the LSL and the UDEC construction of the OLR. Subsequently, the mechanical parameters were assigned separately for the LSL and OLR. For the calculation, the self-weight balance of the LSL in the PFC was calculated firstly. Subsequently, the stress of the bottom wall (contact interface) of the LSL in the PFC was extracted to apply to the top of the OLR in the UDEC to simulate the load of the LSL on the OLR. After balancing the stress calculation of the OLR in the UDEC, the coal seam excavation begins. The excavation distance of each coal seam can be determined based on actual mining. After excavating the coal seam, the mechanical calculations of the OLR in the UDEC was performed. While calculating a certain time step, the top boundary breakage morphology of the OLR was extracted and imported into the PFC to form an LSL bottom wall. This can realize dynamic combining deformation between the OLR and LSL models. In the simulation, a longer calculation time step resulted in higher accuracy but lower efficiency. Mechanical calculations of the LSL were performed based on the new bottom wall in the PFC, and the stress at the bottom wall was extracted to update the stress at the top of the OLR after balancing. The mechanical calculation was continued based on the latest stress at the top of the OLR in the UDEC, and the steps above were repeated until the maximum displacement change at the top of the bedrock was less than 5% (which can be set based on the calculation-accuracy requirements). After calculating the balance, the coal seam was excavated. The specific process is illustrated in Fig. 4.

Combining procedure between two distinct simulation methods

When the contact surface between the OLR and LSL deformed, the breakage morphology of the contact surface was transformed into a “. Dxf” file after certain steps were calculated in the UDEC model. Subsequently, the file was imported into the PFC model using the Fish language such that it became a “wall” for controlling particle deformation. After the PFC model was balanced, the vertical additional stress on the “wall” was extracted and applied to the UDEC model for continuous calculation, which formed a dynamic cycle of combined deformation between the OLR and LSL, thus achieving the dynamic combined deformation simulation of the OLR and LSL.

3.1 Mining and geological conditions of study area

The 12401 longwall face of Shangwan coalmine is the first mining face of 1⁻² Coal Seam, with an advancing length of 5254.8 m, a strike length of 299.2 m, and a design mining height of 8.6 m (maximum of 8.8 m), as shown in Fig. 5a. The overall thickness of the overburden was 120–224 m, and the thickness of the loose layer was 0–27 m. The main roof is fine-grained sandstone with a thickness of 14 m and a distance of 20.2 m from the coal seam. The PKS is a coarse-grained sandstone with a thickness of 22.7 m and a distance of 68.5 m from the coal seam. The strata distribution map is shown in Fig. 5b.

The layout of the longwall face and the occurrence characteristics of rock layers. a 12401 longwall face layout; b The strata distribution

3.2 Model construction and parameter validation

3.2.1 Establishment of OLR model

Based on the geological conditions of the 12401 longwall face of the Shangwan coal mine in the Shendong mining area (Fig. 5b), a two-dimensional discrete-element numerical model of the OLR was constructed, as shown in Fig. 6. The model includes the coal seam floor to the top of the strata, with dimensions of X = 500 m and Y = 228 m. A 100 m coal pillar was implemented at the model boundary to reduce boundary effects. Based on previous DEM simulations (Cheng et al. 2020; Sun et al. 2019; Zhang et al. 2016), rock stratum is composed of rock blocks-joint surfaces-rock blocks. When the overburden was disturbed by mining, the rock blocks can move, rotate, and deform, and nonlinear deformation characteristics such as cracks or slippage can occur on the joint surfaces. Therefore, the rock stratum was discretized into blocks of different sizes. In addition, although Voronoi and Trigon blocks improve the accuracy of calculation for the failure and collapse of intact materials, a greater number of elements improve plastic collapse calculation while simultaneously decrease calculation speed (Le et al. 2018). Therefore, for the entire rock layer, rectangular blocks have higher computational efficiency compared to Voronoi and Trigon blocks (Wu et al. 2019; Gao et al. 2014). Gao et al. (2014; 2022) discovered that a block size of 2.0 m was sufficiently fine to simulate roof caving while maintaining computational efficiency, whereas for thick and hard rock layers, the block size can be increased. Wang et al. (2012) divided the length of a primary key strata block based on the actual weighting step distance. In our study, the average periodic weighting step distance was 18.33 m (Zhang et al. 2021). Thus, the model was divided into blocks with sizes ranging from 2 to 20 m based on the strength of the rock layers. The bottom boundary of the OLR model was fixed in the vertical direction and the top boundary was free. The left and right boundaries were constrained in the horizontal direction. The in-situ stress can be calculated using \(\sigma_{{y}} = \gamma h\) and \(\sigma_{{x}} = K_{{x}} \gamma h\), where K_x is the horizontal lateral pressure coefficient. In this paper, K_x = 1.2 was set based on the measured results of drilling stress relief method.

The model of overburdened strata. a UDEC numerical model; b Mechanical model for numerical analysis of mining site

Additionally, based on the on-site advancing speed, the mining work was advanced 10 m ahead at each excavation step. Five horizontal measuring lines were set in the model, with 50 monitoring points evenly arranged on each horizontal measuring line to record the vertical displacement changes during the mining of the longwall face (as shown in Fig. 6b). In addition, to display the breakage morphology of the bedrock top more intuitively and to further investigate the subsidence law of the OLR, the division of the top block of the OLR is smaller than other rock layer.

3.2.2 Establishment of LSL model

Based on the actual geological conditions of the 12401 longwall face, the dimensions of the LSL model were X = 500 m and Y = 8 m, as shown in Fig. 7. A 100 m coal pillar was implemented at the model boundary to reduce boundary effects. Owing to the significant increase in computational complexity caused by using actual particle sizes in engineering-scale simulations, the particle size must be increased proportionally. For flat-joint or parallel bonding models, owing to the bonding between particles forming a rock mass, the particle size affects the mechanical properties of the rock mass. To eliminate the effect of particle size in the PFC simulation, the ratio of the model size to the average particle diameter should be between 30 and 250 (Cheng and Wong 2020; Liu et al. 2023). However, for the rolling resistance contact model used in this study, the particle size minimally affected the mechanical parameters of each particle (Ai et al. 2011; Scheffler and Coetzee 2023). Therefore, to ensure computational efficiency, the actual particle sizes must be matched as much as possible. Based on a screening test, the particle size range of aeolian sand at the site was 0.3–1 mm. Thus, in this study, to improve the computational efficiency, the particle size was increased correspondingly, with the particle size primarily ranging from 75 to 250 mm. The ratio of the model size to the particle size in this study was larger than those of Park et al. (2005), Liu et al. (2023), and Shi et al. (2021). Based on a drainage experiment, the aeolian sand particle porosity was 33% and the model contained 88,526 particles. To facilitate the observation of particle motion, a trend observation line was set on the surface of the model and measurement points were arranged every 10 m.

The Standard model diagram for numerical simulation

3.2.3 Model-parameter calibration

1. (1)

   Parameter calibration of OLR in UDEC

   The parameters in the block discrete-element simulation primarily included the block and joint parameters. The Mohr–Coulomb plastic model was used in this simulation to determine the mechanical behavior of the block materials, and the Coulomb slip model was used to determine the mechanical behavior of rock contacts between the blocks. The mechanical parameters of the rock layers were determined via experimental laboratory tests. However, the mechanical characteristics between intact rocks on a laboratory scale and rock masses at an engineering scale are different. Therefore, Zhang and Einstein (2004) and Zhang et al. (2019) utilized an RQD index to propose a conversion formula for the mechanical parameters of intact rocks and rock masses.

   $$\frac{{T_{\text{m}} }}{{T_{\text{r}} }} = 10^{0.013\,\text{RQD} - 1.34}$$

   (5)
   $$\frac{{\sigma_{\text{cm}} }}{{\sigma_{\text{cr}} }} = 10^{0.013\,\text{RQD} - 1.34}$$

   (6)
   $$\frac{{E_{\text{m}} }}{{E_{\text{r}} }} = 10^{0.0186\,\text{RQD} - 1.91}$$

   (7)

   where T_r and T_m are the tensile strengths of intact rocks and rock masses, respectively; σ_r and σ_m are the uniaxial compressive strength of intact rocks and rock masses, respectively; and E_r and E_m are the moduli of intact rocks and rock masses, respectively. The mechanical properties of the coal and rock masses after calibration are listed in Table 1. Based on these properties, the laboratory compressive and tensile inversion simulations were performed to determine the joint parameters using a trial-and-error optimization algorithm (Gao and Stead 2014), as listed in Table 2. In addition, the parameter calibration at the engineering scale was conducted, and the deep base-point displacement-meter monitoring data was used to verify the simulation parameters.

   Lithology	Density (kg/m3)	Cohesion (MPa)	Tensile strength (MPa)	Bulk modulus (GPa)	Internal friction angle (°)	Shear modulus (GPa)
   Gritstone	2262	3.9	3.7	7.5	35	3.2
   Sandy mudstone	2254	3.7	3.45	3.4	26	3.4
   Medium sandstone	2275	4.5	2.8	7.7	30	4.7
   Siltstone	2160	5.9	4.5	4.3	25	4.3
   Fine sandstone	2726	4.8	3.9	7.9	33	3.8
   Mudstone	2050	3.3	3.6	3.5	27	2.9
   Coal seam 1–2	1305	2.3	1.53	2.1	24	3.1

   Lithology	Cohesion (MPa)	Internal friction angle (°)	Normal stiffness (GPa)	Tensile strength (MPa)	Shear stiffness (GPa)
   Gritstone	2.36	15	2.3	0.06	1.46
   Sandy mudstone	2.5	12	1.6	0.04	1.36
   Medium sandstone	2.12	13	2.2	0.05	1.32
   Siltstone	2.96	16	2.7	0.02	1.24
   Fine sandstone	2.47	14	2.5	0.06	1.28
   Mudstone	1.4	15	1.2	0.02	0.4
   Coal seam 1–2	1.1	10	1.0	0.01	1.3

2. (2)

   Parameter calibration of LSL in PFC

   The material properties were defined as the inherent characteristics of the particles, such as shape, size distribution, density, Poisson’s ratio, shear modulus, and yield strength. The material-interaction characteristics refer to the characteristics of particles in contact with the boundary surface. Additionally, the rolling and static friction coefficients are important indicators (Jiang et al. 2005; Rorato et al. 2021; Wensrich et al. 2012; Bharadwaj et al. 2010). The friction coefficient between particles determines the repose shape of particles and the static friction coefficient (μ) exerts the most significant effect on the repose angle (Jiang et al. 2005; Goniva et al. 2012; Ai et al. 2011, Coetzee 2020). Therefore, the static friction coefficient between the particles was measured via a slope sliding test, as shown in Fig. 8. The static friction coefficient between the particles was 0.69. However, owing to the small size of windblown sand particles, their rolling friction coefficients could not be measured via mechanical experiments. Therefore, a slump experiment simulation was performed for calibration, as shown in Fig. 9.

   

   Slope sliding test for aeolian sand on the panel

   

   Slump experiment. a Laboratory test; b Numerical simulation

As shown in Fig. 9 and Table 3, grayscale binary processing, hole filling, and boundary extraction were performed on the sand piles from different perspectives. The least-squares method was used to fit the contour boundary curve to a fitting line. The slope of the fitting line is the tangent of the repose angle. Finally, the repose angle of the aeolian sand particles was 26.741°. Meanwhile, as shown in Fig. 10, when the rolling friction coefficient between the particles was 0.47, the simulation test of the solution results showed that the repose angle was 26.31° and the result error was only 1.6%. Thus, based on the simulation, the rolling and static friction coefficients were 0.47 and 0.69, respectively.

Analysis of slump experiment. a Laboratory test; b Numerical simulation results

3.3 Collaborative movement characteristics

3.3.1 Stratified OLR subsidence

Based on the monitoring data of the layered subsidence, the PKS determines the overall deformation of the OLR and LSL. The damage and breakage morphology of the OLR after mining varied with the PKS. The bedrock-top breakage structure is the final manifestation of mining damage during the conduction stage of the OLR under high-strength mining and directly results in surface subsidence. The breakage morphology and subsidence curve during the advancement of the longwall face are shown in Fig. 11.

Bedrock fracture structure and subsidence curve under different advancing distances. a Advancing by 40 m; b Advancing by 100 m

As shown in Fig. 11a, when the longwall face advanced by 40 m, the main roof was initially completely broken and was neatly arranged in the goaf. This is basically consistent with the initial pressure step distance 38 ~ 40 m of the main roof obtained from on-site measurement (Wang et al. 2024). At this point, the maximum subsidence of the main roof was 5074 mm. In addition, the main roof subsidence was relatively large near the open-cut side relative to the longwall face. This is because the contact type between the cut blocks was different from that between the hinged blocks. Specifically, it was face contact, which generated a large friction force under compression to hinder collapse. At this point, the top of the bedrock did not undergo any damage or deformation. As the longwall face continued to advance, the PKS experienced an initial break when it advanced 100 m. As shown in Fig. 11b, the maximum subsidence value of the bedrock top was approximately 2500 mm. Meanwhile, the PKS formed an “three hinge arch” structure due to the bearing capacity of the collapsed roof. After continuing to advance for 15–20 m, owing to their self-weight and compaction effect on the upper strata, the originally stable rock layer rotated continues to undergo rotational fracture, and the broken blocks were arranged neatly in the middle of the goaf due to hinge action, forming a “masonry beam structure”, as shown in Fig. 12.

Full mining breakage morphology

In the simulation, the measuring points were arranged at a distance of 50 m from the open cutting to analyze the dynamic evolution process of the primary control rock layers; additionally, the results obtained were compared with the results measured using the displacement meter, as shown in Fig. 13. The bedrock top underwent severe subsidence until it reached stability at 150 m beyond the measuring point. At this point, the OLR began to enter the initial full-mining stage, and the subsidence coefficients for the main roof and PKS were 0.72 and 0.73, respectively. As it continued to advance and the measuring points entered the full-mining stage, the final subsidence coefficients of the main roof and PKS were 0.74 and 0.75, respectively. Because the longwall face was located in a thick overburden area, under the action of the overlying strata load, the dynamic subsidence heights of the main roof and PKS were collaborative. When the full-mining stage was reached, owing to the effect of the coefficient of fragmentation and expansion in the caving zone, the subsidence coefficients of displacement meters 9 (main roof) and 5 (PKS) were 0.75 and 0.73, respectively. The dynamic subsidence trend of the primary control rock layer in the numerical model was the same as the on-site measurement trend, with an error of only 1.3%–2.7%, thus verifying the reliability of the OLR numerical model.

Comparison of dynamic settlement points of main control rock layers. a Primary key stratum; b Main roof

The subsidence characteristics indicate the decisive role of the PKS. Based on the occurrence characteristics, although the PKS controls the OLR breakage morphology, a thick and soft buffer layer measuring approximately 40 m was indicated between the PKS and bedrock top. After the PKS was broken, upward damage conduction decelerated. Meanwhile, the movement characteristics of the bedrock top were the ultimate external manifestations of mining damage as well as the “bridge” that conducted mining damage to the LSL, thus directly causing dynamic movement and distribution under the effect of mining. Therefore, the movement law must be analyzed. Additionally, a measuring point P was selected at a distance of 100 m from the open-off cut at the top of the model to analyze the dynamic subsidence process, as shown in Fig. 14.

Subsidence process of bedrock top. a Numerical simulation results; b Evolution law of separation layer; c Verification of simulation results

As shown in Fig. 14a, when the longwall face advanced 100 m, the surface of the bedrock top initially presented an “asymmetric inverted V-shaped” subsidence trough. When it advanced 150 m, the surface subsidence trough reached its maximum depth. At this time, the surface entered the initial full-mining stage. As it continued to advance, the central region of the goaf become more compact and the subsidence trough appeared flat gradually. When it advanced 240 m, the surface subsidence trough developed completely and expanded only laterally. At this time, the surface entered the full-mining stage.

Based on the monitoring of the strata subsidence system, when the longwall face entered the initial full-mining stage, a separation layer measuring approximately 1444 mm appeared at the overburden. In addition, because of the compaction and self-recovery of the overburden fracture, the bedrock remained in a slower compaction state until the internal separation layer was closed, and no movement space was available for a long time after it entered the stable stage (Cai et al. 2023; Hu et al. 2022). Based on the actual process of a fully balanced numerical simulation and the short period of on-site measurement, the separation layer and the existing subsidence amount at the bedrock top were combined and considered as the final ultrafull mining subsidence result, as shown in Fig. 14b.

When it advanced through the measuring point by approximately 140 m, the OLR was in the active movement stage, and the numerical model during this stage was consistent with the movement trend of displacement meter 1 obtained from field measurements. When it entered the fully compacted state, the maximum subsidence coefficient was 0.77, whereas that of the numerical model was 0.74, with an error of only 3.9%. This indicates the reliability of the numerical model and the rationality of the simulation method for modeling the effect of the bedrock top-breakage morphology on the LSL movement state.

3.3.2 Characteristics of collaborative surface movements

Based on the simulation plan shown in Fig. 4, the bedrock top displacement (Fig. 14) was imported into the PFC^2D simulation software for secondary stability calculation, and the results are shown in Fig. 15.

Surface subsidence cloud map and morphology

As shown in Fig. 15, “step-like” subsidence and collapse pits appeared on the surface of the open cut and stope areas, with a step drop of 1.5–2.5 m and a collapse pit depth of approximately 1 m. In addition, based on the subsidence cloud map, the subsidence basin of the OLR-LSL combination presented central symmetry. The maximum subsidence value in the middle of the goaf was 7–8 m. In the LSL model, a few particles sank to a height of 9 m. This is because aeolian sand particles possess strong fluidity, which results in their irregular flow during subsidence and individual particles subsiding at locations higher than the coal seam. The surface-measurement point data of the model were retrieved and processed, and the surface-subsidence curve is shown in Fig. 16.

Surface and bedrock top subsidence curve. a Overlying rock and surface subsidence curve; b Comparison diagram of surface subsidence

As shown in Fig. 16a, by comparing the bedrock-top movement data with the surface-subsidence values, it was discovered that, owing to the complexity of the particle movement, the data at each measurement point were relatively chaotic. However, they were consistent with the subsidence of the bedrock top and presented an “inverted trapezoidal” subsidence trough. Compared with the bedrock top subsidence, the surface subsidence decreased, thus indicating that the LSL imposed a certain absorption effect on mining damage. The maximum subsidence point of the bedrock top was 70 m from the open-off cut, with a maximum subsidence value of 7207 mm and a subsidence coefficient of 0.8. Meanwhile, the maximum subsidence point of the overburden-loose layer combination was 75 m from the open-off cut, with a subsidence value of approximately 6637 mm and a subsidence coefficient of 0.74, which corresponded to decrease by 8%. This indicates that the absorption of mining damage by the LSL is quantitatively manifested as a decrease in surface subsidence compared with the bedrock-top subsidence value.

Based on the monitoring results of surface-displacement observation lines (see Fig. 16b), during the ultrafull mining stage, the results yielded by the “subsidence trough flat bottom” and LSL models were consistent, whereas the results yielded by the numerical simulation method on both sides of the subsidence trough were relatively slow, although the overall trend was consistent.

4.1 Field monitoring methods

To effectively study the characteristics of overburden-surface subsidence in the extra-high mining face of Shangwan Coalmine, and grasp the conduction law of mining damage, the OLR-LSL chain damage subsidence monitoring method was constructed. For the control and production of the PKS of the OLR, a strata subsidence monitoring system was installed in the SD₁ pre-mining hole, and a total of 9 displacement meters were installed at 126, 110, 99, 88, 71, 52, 43, 34, and 26 m away from the coal seam, as shown in Fig. 17a. Meanwhile, two displacement observation lines (K line and L line) are arranged at a distance of 300 m from the open-off cut along the centerline of the working face on the surface, with a total length of 1000 m and a spacing of 20 m between control points, as shown in Fig. 17b.

On-site monitoring layout. a Displacement meter arrangement; b Layout of surface movement observation stations

4.2 Monitoring results and analyses

1. (1)

   Overburden strata subsidence and damage law

   The process of advancing the longwall face 300 m past the SD₁ hole is used as the analysis range for the conduction characteristics of mining damage in the overburden. The monitoring results of the strata subsidence system reflect the dynamic movement law of thick overburden during the conduction process of mining damage. The relative surface subsidence value curve of each measuring point is shown in Fig. 18.

   

   Subsidence curve of OLR relative to the surface

   As shown in Fig. 18, during the active stage, the subsidence curves of the main roof and its control strata exhibit a “step-like” subsidence feature. When the longwall face advances through the borehole for 12 m (advancing 20 m), displacement meters 5–9 suddenly experience a slight disturbance, and the relative surface subsidence only suddenly increases to 280 mm. When it advances 43 m past the borehole (advancing 51 m), there is a large-scale collaborative subsidence of the displacement meters 3–9 relative to the surface. At this time, the PKS was completely broken. When continuing to advance through the borehole for 100 m (advancing 108 m), displacement meters 6–9 experiences severe settlement, with a maximum relative surface subsidence of 1600 mm and displacement meters 3–5 only subsided 1000 mm. Especially when advance through the borehole for about 70 ~ 80 m, the displacement meter 9 (main roof) suddenly subsided by 900 mm. At this time, displacement meters 3–5 only subsided by about 150 mm, which indicates that the PKS subsidence trend was relatively slower. From this, it can be inferred that the main roof and its control strata have undergone collaborative periodic cutting, while the PKS has also undergone bending.

   Due to the collaborative subsidence of the accompanying strata above, it loaded to compress the collapsed blocks in the goaf, resulting in collaborative large-scale subsidence of the displacement meters 3–9. When it advances through the drilling hole for 140 m (advancing 148 m), the OLR reaches initial-full subsidence, and the subsidence value at this time is as high as 1064 mm (displacement meter-4) ~ 1618 mm (displacement meter-7); It is worth noting that the relative surface subsidence of displacement meters 1 and 2 is only 75–126 mm, which exists a separation space with a height of approximately 1131 mm (Fig. 14b). Therefore, they did not collaborative subsidence with the breaking of the PKS. Therefore, based on the geological conditions (Fig. 3b), the borehole opening (surface), displacement meter-1, and displacement meter-2 are jointly controlled by the PKS and the thick-soft rock buffer layer.

2. (2)

   Surface subsidence and damage law

   During the process of advancing, continuous observations were made on the surface K-line (strike line) and L-line (dip line), with the K-line subsidence curve shown in Fig. 19a. When the longwall face advancing 100 m, a subsidence trough appears on the surface. As the longwall face continued to advance, the surface movement observation area experienced strong subsidence and deformation when it advanced by about 350 m. The subsidence basin was developed and formed, forming an “asymmetric V-shaped” shape. At this time, the maximum subsidence value (K₁₉) had reached 6315 mm, with a maximum subsidence coefficient of 0.72, as shown in Fig. 19b. During the continuous advancement process, the subsidence basin develops horizontally as an “asymmetric inverted trapezoid”. It is consistent with the OLR-LSL combining numerical simulation results (Fig. 16b). When it advances about 700 m, the surface subsidence basin is completely stable, and the surface enters the “ultra-full mining stage”.

   

   Observation of surface movement. a Surface subsidence curve of K line; b K₁₉ dynamic subsidence curve

Surface deformation is the result of combining the characteristics of fracture morphology in the OLR with the response motion of the LSL. Previous studies showed that the DEM can effectively simulate the physical deformation of rocks and granular materials, and thus effectively simulate the large deformation characteristics of mining activities (Chen and Hu 2018; Sun et al. 2019; Wang 2018; Khanal et al. 2017). Based on the differences in the physical and mechanical properties and the nonlinear deformation characteristics between OLR and LSL materials, a combined simulation method of block + particle discrete elements was proposed to address the limitations of the UDEC and PFC software in simulating OLR and LSL materials. This combining simulation method achieves a dynamic combining deformation simulation between the OLR and LSL through stress transfer in particle models and via the cyclic conduction of breakage morphology at the contact surface between the OLR and LSL. The results were compared and analyzed with the results of on-site bedrock deep base-point displacement meter monitoring and surface-subsidence monitoring, and the reliability of the results was verified.

Compared with individual UDEC simulations, although rigid blocks can effectively simulate the fracture of OLR and the development and evolution characteristics of mining-induced fractures, they can only slide along the predetermined fracture surface (Gao et al. 2014), which is not suitable for simulating the movement and deformation characteristics of LSL particles. Compared with individual PFC simulations, although the particle model can reproduce the development characteristics of the OLR fractures well, the simulation effect of the rock fracture structure is unsatisfactory. In addition, because of the limitation of the model size, if the model is fully modeled under a certain particle size, then a lower computational efficiency will be resulted (Zhao and Wang2011; Feng and Owen 2014).

However, in our simulation method, the block size in the UDEC directly affects the fracture morphology of the OLR. Thus, the block size should be finely divided for key areas to ensure computational accuracy; however, this reduces the computational efficiency. The same issue applies to the particle size in the PFC simulations. In addition, the moisture content, particle geometry of the loose layer with different burial depth are different, it significantly affects the particle parameters (Zhou et al. 2022; Sadek et al. 2011). Therefore, further research is required to simulate the LSL under different burial depth. Moreover, our combining simulation is indirectly implemented by the cyclic conduction of the breakage morphology of the contact surface and the additional vertical stress of the LSL in certain calculation steps. Hence, the interval time step directly affects the simulation accuracy. The smaller the number of interval time steps, the better the simulation accuracy. But this will significantly increase the computational consumption. In this paper, the appropriate number of interval time steps were determined by comparing with on-site measurement results. In addition, the LSL thickness simulated in this study is relatively small, and causes based on thick LSLs require further verification. Thus, more engineering case studies shall be conducted in the future.

In this paper, a novel block–particle discrete-element simulation method was developed that matched the double medium of OLR and LSL in coal mining. And two discrete element numerical simulation methods, UDEC and PFC2D, were used to construct models for OLR and LSL. By combining the stress of the LSL to the loading of the OLR and the breakage morphology of the bedrock top to control the LSL movement, the mining collaborative damage reproduction of the double medium was achieved. Meanwhile, combined with on-site measurements, the characteristics of fracture and subsidence of the bedrock and surface were revealed. The main conclusions are as follows:

1. (1)

   A linear rolling-resistance contact model was applied to loose-layer particles to match the mutual hindrance between irregularly shaped particles of aeolian sand. Static and rolling friction coefficients were calibrated via the slope sliding method, a slump experiment, and trial-and-error simulations. The static and rolling friction coefficients were determined to be 0.69 and 0.47, respectively.

2. (2)

   The combining simulation method was used to elucidate the effects of the bedrock-top breakage morphology and LSL physical characteristics on surface subsidence. The surface subsidence basin was consistent with the bedrock-top breakage morphology and presented an “asymmetric inverted trapezoid”. The maximum subsidence coefficients of the surface and bedrock top were 0.74 and 0.8, respectively, which explains the absorption effect of the LSL on mining damage. Additionally, the surface formed “step-like” subsidence and collapse pits at both boundaries of the subsidence basin.

3. (3)

   The OLR was guided by the “double-control layer and thick-soft rock buffer layer” and indicated “cut-off grouping” subsidence, whereas the surface formed collaborative subsidence with the thick and soft rock buffer layer. After reaching the full mining stage, the subsidence trough extended horizontally in an “asymmetric inverted trapezoid” shape as the operating face advanced. The maximum surface subsidence value (K₁₉) was 6315 mm, and the maximum subsidence coefficient was 0.72. Meanwhile, the maximum subsidence coefficient of the bedrock top was 0.77, with relative errors of 2.8% and 3.9% in the numerical simulation compared with the field measurement. Thus, the reliability of the block–particle discrete-element combining simulation method was verified.