The Kunyuan Coal Mine belongs to Daliuta Town, Yulin City, Shaanxi Province, China, and is located in the northern part of the Loess Plateau in Northern Shaanxi Province, at the southeastern edge of the Maowusu Desert. The geomorphological unit is a typical loess hilly area, and the surface is basically covered by wind-bound sand. The gullies are deep in the bedrock, and the branch gullies are very well developed. There are two gullies crossing the surface of 113101 working face from northwest to southeast (Figure 1). Gully No. 1 is a wash gully, with average slope angles of 35.82° and 35.30° for the downhill and uphill sections, respectively (the downhill section is the side of the slope where the thickness of the overburden rock gradually decreases along the direction of advancement, and the opposite for the uphill section). There is a large weathering zone and burned rock on the east side of the slope and water flow in the gully all year, with an average flow rate of 197.30 m³/h. Gully No. 2 is a cut gully, with average slope angles of 28.44° and 21.59° for the downhill and uphill sections, respectively. There is water flow only in the rainy season. Long-term water erosion has destroyed the structural integrity of the overburden rock, and two weathering erosion zones with vertical depths of 0.2–4.3 m have formed along the riverbed.

The 113101 working face is mainly mining a 3^–1 coal seam with a buried depth of 29.36–91.80 m. The lowest depth of the coal seam in Gully No. 1 is 7.49 m, compared to 16.42 m in Gully No. 2. The average thickness of coal seam is 1.81 m, the structure is simple, and the bottom contains a stratum of gangue. The immediate roof is gray, gray–white mudstone, and sandy mudstone, with mud cementation and lamination development, containing a small amount of mica and plant fossils, and belongs to class III immediate roof; the main roof is gray and gray–white dense coarse-grained sandstone, with mud cementation and wavy lamination development, containing plant fossils, and belongs to class II main roof; and the floor is mainly gray and gray–white dense coarse-grained sandstone with mud cementation. The physical and mechanical parameters of the coal rock are shown in Table 1. The working face is 776 m in trend length and 280 m in inclination length. A longwall comprehensive mechanized coal mining method and natural collapse method are used to manage the roof, with a mining height of 1.80 m, five cycles of daily coal cutting, and a cut-off depth of 0.8 m. The purpose of coal water co-mining at working face 113101 is to recycle and utilize mine water, divert surface rivers, and restore the surface ecology of the mine area as much as possible after mining. As shown in Figure 2, the exploitation of water resources is divided into two parts: damming the surface to intercept the flow and pumping water in the collecting drain of goaf. The dam is located upstream of the first gully, 50 m from the tailgate of 113101 working face. The water collecting drain is advanced in the floor of the goaf in the middle of Gullies No. 1 and 2 and the water collecting drain is advanced in the floor at 30 cm in height and covered with a six-mesh screen.

Taking the Kunyuan Coal Mine 113101 working face as the engineering background, according to the engineering, geological, and production conditions as well as the physical and mechanical parameters of coal rocks (Table 1), a 3D model of the surface gully was established by Rhino software. Then, the grid was divided by Griddle and imported into FLAC3D software to generate a 3D numerical simulation model (Figure 3). The model size was L × W × H = 620 m × 580 m × 102 m, with 144,853 grid points and 795,173 cells. The Mohr–Coulomb principal structure model with large deformation stress mode was used. Taking the B–B section as reference, we analyzed the overburden stress evolution characteristics of shallow buried coal seam mining through the downhill section of Gully No. 2 (X = 190–270 m), uphill section of Gully No. 2 (X = 270–330 m), downhill section of Gully No. 1 (X = 330–390 m), and uphill section of Gully No. 1 (X = 390–535 m) to reveal the overburden transport law.

The maximum principal stress distribution of the B–B section when the working face advances to Gully No. 2 (X = 190–332 m) is shown in Figure 4. The direction and size of the maximum principal stress deflection indicates that the overburden rock above the mining area formed a stress arch structure, with the rock strata below the stress arch in a stress release state and the rock strata above the “stress arch” in a compressive stress concentration state. This stress distribution state causes the surface behind the working face to be subjected to tensile stress, which is expressed as extrusion, and the surface above the working face to be subjected to tensile stress, which is expressed as tension. Specially, the bottom of Gully No. 2 was subjected to strong tensile stress effect for a long time and the effect was obvious, and the details are as follows.

The stress transfer law in Gully No. 1 was approximately the same as that in Gully No. 2. The details of the dynamic changes of the maximum principal stress during the advance of Gully No. 2 were as follows. As shown in Figure 5, the working face was advanced to X = 190 m. According to the route of maximum principal stress transmission, the stress arch was in the shape of a “saddle,” similar to the shape of the surface, the overburden rocks above the stress arch were subject to compressive stress, and the overburden rocks below the stress arch were subject to tensile stress. From the surface stress situation, the back of the working face was subjected to compressive stress and the front of the working face was subjected to tensile stress, in which the concentration of tensile stress at the bottom of Gully No. 2 was the most obvious, with a maximum value of 1.83 MPa.

As shown in Figure 6, when the working face was advanced to 272 m, the stress arch shape returned to normal. Taking the gully bottom as the boundary, the surface behind the working face showed obvious compressive stress concentration and that in front of the working face showed tensile stress concentration. Currently, the most complicated changes were in the gully bottom and the magnitude of stress changes was large, making it easy for the overburden rock to collapse with obvious elevation difference in the gully bottom.

In summary, the macroscopic stress transfer path was as follows: first, the stress behind the working face was transferred to both sides along the trend and the stress transferred to the working face was arched, resulting in the concentration of compressive stress behind the working face; then, the stress in front of the working face was transferred to the working face, and the tensile stress was concentrated in front of the working face; finally, the stress in the downhill and uphill sections above the goaf was transferred to the gully bottom.

When the working face was advanced to X = 280 m, the horizontal extrusion stress was concentrated in the lower part of the downhill section of Gully No. 2 and the gully bottom, the maximum horizontal extrusion stress was 4.68 MPa, the maximum horizontal tensile stress in the uphill section was 1.53 MPa. The downhill section of Gully No. 2 was above the goaf, and sliding fractures would be formed throughout the downhill section while extrusion fractures would only be formed at the base of it. The uphill section was above the coal and contained over-extension fractures as shown in Figure 7A. When the working face was advanced to X = 340 m, the horizontal extrusion stress at the bottom of Gully No. 2 was 8.00–12.00 MPa and the stress state changed from tensile stress to minor stress to compressive stress. The fractures healed or changed types to form other fractures, as shown in Figure 7B. When the working face was advanced to X = 410 m, a small range of horizontal compressive stress was concentrated in the lower part of downhill section of Gully No. 1, the maximum compressive stress was 6.00 MPa, the formed fractures were easily healed to varying degrees, a small range of horizontal compressive stress was concentrated in the lower part of the uphill section of Gully No. 1, the maximum horizontal tensile stress was 1.75 MPa, and tensile fractures were prone to form along the weak surface of rock joints, as shown in Figure 7C. When the working surface was advanced to X=460 m, horizontal extrusion stress concentration occurred at the bottom of Gully No. 1, and extrusion and collapse fractures formed at the base of the slope. The loose layer and weathering zone of the uphill slope body of Gully No. 1 were staggered with the bedrock under the action of tensile stress, which made it easier for sliding fractures to form, as shown in Figure 7D.

The horizontal stress of the surface first increased and then decreased during across-gully mining. The horizontal stress of the surface was 2.11–4.02 MPa in the downhill section and 5.92–17.35 MPa in the uphill section. The surface of the gully was in a horizontal compressive stress state far from the working face, as shown in Figure 8A. The vertical tensile stress of the surface first increased and then decreased. The vertical tensile stress at the surface of Gully No. 2 was 0.69–3.74 MPa, with the vertical stress at the surface being the largest across the gully. The vertical tensile stress at the surface of Gully No. 1 was 0.69–1.56 MPa and the uphill section was in a vertical compressive stress state behind the working face, as shown in Figure 8B. The vertical tensile stress reflected to the overburden collapse state and the tensile stress range in front of the working face included the overburden fracture development and overbreak areas. The tensile stress range in the gully area included the collapse fracture development area. The bottom of Gullies No. 2 and 1 was in a high-tensile stress state, which encouraged collapse fractures.

To verify the numerical simulation results and further determine the dynamic development process of fractures in across-gully mining, a 2D physically similar simulation model was established based on the B–B section (Figure 9), with a geometric similarity ratio of 1:200, capacity similarity ratio of 0.64, stress similarity ratio of 3.2 × 10^–3, and similar material ratios for each rock strata, as shown in Table 2. The model dimensions were L × W × H = 2000 × 200 × 506 mm. The coordinate system of the physically similar simulation model and the numerical simulation model were the same, and the coordinates of the lower left corner of the physically similar simulation model were X = 140 m and Z = 20 m) The left end of the model had no contact with the model frame and the model was advanced 20 cm to simulate the goaf in front of gully to avoid the error caused by the left boundary of the model.

When the working face advanced to X = 230 m, tensile fracture F2 developed to the top downhill section of Gully No. 2, the overburden breakage angle was 120°, and the topmost loess of the model was destroyed into seven parts of different sizes, with length ranging from 3.0 to 5.6 cm (Figure 10A). When the working face advanced to X = 250 m, the immediate roof collapse and main roof breakage caused the overburden rock to form a big intact rock. The center of gravity of this rock mass shifted above the mining areas, which led to the expansion of tensile fracture F4, the formation of tensile fracture F3, and the narrowing of fracture F2 (Figure 10B); fracture F4 was located at the base of the downhill section of Gully No. 2, which connected Gully No. 2 with the goaf and became a water–sand inrush channel. The formation of the big intact rock indicated that the overburden load decreased after the overburden thickness decreased, which led to the change of cycle collapse distance, and fractures were formed in the thin bedrock with low strength and weak surface as the working face advanced.

As shown in Figure 11, when advancing to X = 270 m, the working face was located below the bottom of Gully No. 2. The surface rift F4 broke downward, forming a break line of approximately 5 cm long, and stopped developing at 11.2 cm above the working face. In addition, a new break line developed from the bottom of Gully No. 2 downward, which exceeded the working face by 10.3 cm, called an overburden overbreak line. This line stopped developing at 6.1 cm above the working face and was the prototype of fracture F5. When the working face advanced to X = 280 m, the overburden overbreak line developed into collapse fracture F5, which penetrated the overburden to the working face from the surface. At this point, the overburden rock was cut off above or in front of the hydraulic supports, which caused the collapsed roof to press on and damage the supports. The width of fracture was approximately 1.8 mm, the elevation difference between both sides of fracture was approximately 2.6 mm, and a fracture with width of 2.3 mm and elevation difference of 2.6 mm was formed on the surface. In addition, the secondary development of fracture F3 occurred with the formation of fracture F5. When the working face advanced to X=290 m, tensile fracture F6 was formed and the width and elevation difference of fracture F5 were reduced to 0.6 and 1.1 mm, respectively. Fracture F6 was formed by tensile damage but showed the characteristics of collapse fracture at its surface. The width of surface fracture was 1.5 mm and elevation difference was 0.9 mm, which was presumed to be caused by the rotation of broken overburden to goaf. When the working face was advanced to X = 310 m, tensile fracture F7 was formed, and fractures F5 and F6 and the surface fractures were basically healed.

As shown in Figure 12, when the working face advanced to X = 330 m, it entered the downhill section of Gully No. 1 and tensile fracture F8 was formed, forming fractures on the surface. When the working face advanced to X = 350 m, the upper part of fracture F8 narrowed, the surface fractures healed, and the lower part showed secondary expansion; tensile fracture F9 was formed in front of the working face, which overshot the working face by 6.4 cm and was located at the base of the downhill section of Gully No. 1; another big intact rock was formed between fractures F8 and F9. The load of this rock was transferred to supports, which crushed the supports. When the working face advanced to X = 370 m, fractures F8 and F9 narrowed and collapse fracture F10 was formed at the base of the downhill section, which went straight from the surface to the working face and caused support crushing and water–sand inrush. In addition, the fracture F10 formed a surface crack with an elevation difference of 0.37 mm. The left side of this surface crack sank, and the right side was elevated, that is, the left side showed a collapse fracture and the right side showed an extrusion fracture. It is presumed that the overburden rock was caved and formed a collapse fracture, then the broken rock strata and the big intact rock in the back goaf deflected toward the working face, resulting in the left side of fracture F10 squeezing the right side and creating a collapse-extrusion surface crack on the surface near collapse fracture. When the working face was advanced to X = 390 m, tensile fracture F11 was formed, the overburden elevation difference between the two sides of fracture F10 was obviously reduced, and the surface fracture was healed to some extent; however, the extruded and lifted rock on the right side was still higher than the surrounding surface.

As shown in Figure 13, when the working face advanced to X = 410 m, tensile fracture F12 was formed and fracture F11 showed no obvious change. When the working face advanced to X = 430 m, tensile fracture F13 was formed, fracture F11 was narrowed, fracture F12 mostly healed, and both surface cracks disappeared. When the working face advanced to X = 450 m, tensile fracture F14 was formed, fracture F13 shrank, and the surface fracture healed. When the working face advanced to X = 470 m, a large collapse space was formed after the main roof collapsed, which led to further collapse and compaction of the collapsed rock strata on the left side of fracture F14. Fracture F14 was transformed into a collapse fracture and the elevation difference between the two sides of the fracture was 0.22 mm, forming a tensile-collapse crack on the surface.

The fractures in different areas had distinct characteristics, shown in Table 3. The fractures healed more easily in the gully bottom, while the junction between the base downhill section and gully bottom was most likely to collapse; the fracture formation of the slope body with obvious change of overburden thickness tended to be lagging behind and the fracture formation in the downhill section tended to form a big intact rock after the fracture formation was lagging behind. The fracture formation lagging behind in the uphill section tended to trigger the secondary development of the existing fractures and turned into collapse fractures. The distribution of each type fracture and the dynamic development of representative fractures are shown in Figure 14.

Possible disasters that can occur during the mining of the 113101 working face are water–sand inrush, roof caving and supports crushing, and gully-induced landslide. The formation mechanisms of these disasters are as follows.

The disaster areas of the 113101 working face crossing gully was divided according to the types of disaster-causing fractures and disaster-causing mechanism as follows: water–sand inrush areas, water and sand covered areas which easily form water–sand inrush channels, roof caving and support-crushing areas like the gully bottom and the base of the slope where collapse fractures easily formed, landslide I areas where the displacement of the middle and upper part of the slope section evolves slowly and landslides easily occur along the weak surface after a large area of hanging roof, and landslide II areas where the inverse slope section is prone to slope toe instability landslides. These disaster areas are shown in Figure 17.

The disasters caused by gully mining are interrelated. Therefore, they should be analyzed systematically and comprehensively for integrated disaster prevention and control, and corresponding disaster prevention and control measures should be formulated, as shown in Figure 18. The following measures can be taken for the effective prevention and control of disasters:

(1) Water–sand inrush prevention measures: for all types of fractures, the core of prevention is to remove water and sand and close the water–sand inrush channel.

In addition to the aforementioned targeted prevention and control measures, other methods can assist in the prevention and control of various disasters, such as strengthening drainage, advancing faster with small cutting depth, reducing mining height, and moving supports with pressure. In the production process, monitoring means and prevention awareness should be strengthened. The support of the working face and roadway should be strengthened moderately. The mining scheme can be further optimized and a comprehensive command, disaster prevention, and control system can be established.

During the across-gully mining of the 113101 working face in the Kunyuan Coal Mine, the aforementioned prevention and control measures were taken. No disasters occurred during the across-gully mining task, which ensured the safety of personnel and normal production and demonstrated the reliability of these measures.