Study on Overlying Strata Movement and Surface Subsidence of Coal Workfaces with Karst Aquifer Water

: The overlying strata layers of coal workfaces with karst aquifer water normally causes serious safety problems due to the precipitation, drainage and water inrush, such as a wide range and long term of surface subsidence. In this study, by taking 10,301 working faces of the Daojiao coal mine in Guizhou Province as the engineering background, the numerical model of water-bearing strata with ﬂuid-solid coupling was established by using UDEC to illustrate the laws of overlying strata movement and surface subsidence. A theory model was proposed to calculate the surface settlement caused by the drainage of aquifer based on the principle of effective stress modiﬁed by the Biot coefﬁcient α b . The results showed that the corresponding maximum value (0.72 m) and the range of the surface subsidence with the occurrence of karst aquifer water were larger than that of the overlying strata without karst aquifer water (e.g., the maximum value of surface subsidence with 0.1 m). Moreover, the surface subsidence caused by the drainage of aquifer accounted for 17.8% of the total surface subsidence caused by coal mining. According to the ﬁeld monitoring of surface subsidence in 10,301 working faces, the maximum value was 0.74 m, which was highly consistent with the results of numerical simulation and theoretical analysis. It veriﬁed the accuracy and reliability of the numerical model and the theory model in this study.


Introduction
As one type of geological disaster, surface subsidence generally occurs slowly and is difficult to detect, and once the surface has subsided it is almost impossible to completely recover. The ecological environment has obviously deteriorated due to surface subsidence. It has caused serious safety problems and economic losses in industrial and agricultural production, transportation and people's lives.
The whole dynamic movement, deformation and failure process of overlying strata in the period of coal mining was complex due to the special conditions of coal bed and the particularity mechanism of coal seam failure. Therefore, the accuracy regarding the overlying strata movement and surface subsidence played a vital role to ensure coal mining with safety, efficiency, intelligence and green practices. Normally, the overlying strata can be divided into three different moving zones in longwall mining, and the general understanding of the three zones in the horizontal and vertical directions have been proposed, respectively [1][2][3][4][5][6]. Liang et al. established the ANSYS numerical model to illustrate the surface-subsidence laws of thick alluvium caused by coal mining [7]. Xu et al. revealed that the number of key overburdened strata in deep mining was generally larger than that in shallow mining, and that the distance between the main key overburdened strata and the mining coal seam in deep mining was generally greater than that in shallow

Engineering Background
The Daojiao coal mine is located in Songkan Town, Tongzi County, Guizhou Province with geographical coordinates 106 • 53 49 -106 • 54 25 E and 28 • 30 14 -28 • 31 32 E. The designed production capacity of the technical transformation of the mine was 300,000 t/A. C 3 as the main coal seam in the Daojiao coal mine was located in the middle and upper part of Longtan Formation, about 26 m away from the limestone of Maokou Formation, and the thickness of the coal seam was in the range of 1.86 to 2.26 m with an average thickness of 2.02 m. In addition, the dip angle and the average buried depth were 7 • and 210 m, respectively. The 10,301 working faces were laid out in this main coal seam, and the length, mining height and advancing length were 160 m, 2.02 m and 540 m, respectively. Moreover, the roof lithology in the 10,301 working faces was mudstone, carbonaceous mudstone and locally intercalated siltstone, while the floor lithology was clay rock and mudstone as shown in Figure 1. The main aquifers and impermeable stratum in the mining area were Quaternary (Q) pore aquifer, impermeable stratum of Triassic system (T 1 y 3 , T 1 y 1 ), aquifers of Triassic system (T 1 y 2 ), impermeable stratum of Permian (P 3 l) and aquifers of Permian (P 3 c, P 3 c).

Engineering Background
The Daojiao coal mine is located in Songkan Town, Tongzi County, Guizhou Province with geographical coordinates 106°53′49″-106°54′25″ E and 28°30″14″-28°31′32″ E. The designed production capacity of the technical transformation of the mine was 300,000 t/A. C3 as the main coal seam in the Daojiao coal mine was located in the middle and upper part of Longtan Formation, about 26 m away from the limestone of Maokou Formation, and the thickness of the coal seam was in the range of 1.86 to 2.26 m with an average thickness of 2.02 m. In addition, the dip angle and the average buried depth were 7° and 210 m, respectively. The 10,301 working faces were laid out in this main coal seam, and the length, mining height and advancing length were 160 m, 2.02 m and 540 m, respectively. Moreover, the roof lithology in the 10,301 working faces was mudstone, carbonaceous mudstone and locally intercalated siltstone, while the floor lithology was clay rock and mudstone as shown in Figure 1. The main aquifers and impermeable stratum in the mining area were Quaternary (Q) pore aquifer, impermeable stratum of Triassic system (T1y 3 , T1y 1 ), aquifers of Triassic system (T1y 2 ), impermeable stratum of Permian (P3l) and aquifers of Permian (P3c, P3c).

geological characteristics
Purple red, yellow clay and gravel, 0-7 m thick Purple red, dark purple with yellow green mudstone, with thick layered marl, limestone. Thickness greater than 50 m Deep gray, gray thin to medium thick interlayers micro to fine crystalline limestone, argillaceous limestone and argillaceous limestone, the bottom is yellow green, brown gray shale, calcareous s h a l e , c a l c a r e o u s m u d s t o n e . 1 2 0 ~ 1 4 0 m t h i c k .
Mudstone and argillaceous limestone, the bottom is yellow green, brown gray shale, calcareous shale, calcareous mudstone.
The lithology is gray-deep gray medium-thick layer to fine-grained gray. Bioclastic limestone, containing a small amount of black carbon mudstone, bottom with 0.10 ~ 0.50 m black mudstone. The thickness is about 50m. Integrated contact with underlying strata.
The lithology is mainly composed of gray, dark gray mudstone, silty mudstone, clay rock, pyrite-bearing clay rock and coal, with thin lenticular siderite ( limonite after surface weathering ). Carpiopods, plant fossils. The thickness is about 60m. Contact with underlying strata false integration.
Light gray, dark gray medium to thick layered limestone. Bioclastic limestone, lenticular rock in the middle of the gray, dark gray thick to full thick block fine limestone. The upper part is light gray, gray thick massive limestone. Thickness is greater than 100m.

Model Establishment and Parameters Selection
By taking 10,301 working faces of the Daojiao coal mine in Guizhou province as the engineering background, numerical simulation models using fluid-solid coupling with the length and height of 500 m × 240 m, respectively, were established to illustrate the influence of coal mining on the overlying strata movement and surface subsidence while both considering karst aquifer and not, as shown in Figure 2. In this study, the Mohr-Coulomb theory was adopted as the failure yield criterion and the contact Coulomb slip model was selected for the structural plane. Moreover, the upper part of the model was the free boundary and the lower part of the model was restrained vertically with lateral constraints on the left and right boundaries. In terms of overlying strata containing two extremely thick karst aquifers, a smaller block in these layers was set because its integrity and strength were porous. In addition, there was an impermeable stratum of mudstone between the two karst aquifers with good integrity and high strength to set as the larger block, and the block size of other overlying strata gradually increased with the distance away from coal seam. Moreover, the 50 m coal pillars were reserved on the left and right to eliminate the boundary effect with the working face total advancing 400 m.

Model Establishment and Parameters Selection
By taking 10,301 working faces of the Daojiao coal mine in Guizhou province as the engineering background, numerical simulation models using fluid-solid coupling with the length and height of 500 m × 240 m, respectively, were established to illustrate the influence of coal mining on the overlying strata movement and surface subsidence while both considering karst aquifer and not, as shown in Figure 2. In this study, the Mohr-Coulomb theory was adopted as the failure yield criterion and the contact Coulomb slip model was selected for the structural plane. Moreover, the upper part of the model was the free boundary and the lower part of the model was restrained vertically with lateral constraints on the left and right boundaries. In terms of overlying strata containing two extremely thick karst aquifers, a smaller block in these layers was set because its integrity and strength were porous. In addition, there was an impermeable stratum of mudstone between the two karst aquifers with good integrity and high strength to set as the larger block, and the block size of other overlying strata gradually increased with the distance away from coal seam. Moreover, the 50 m coal pillars were reserved on the left and right to eliminate the boundary effect with the working face total advancing 400 m. According to the field-monitoring hydrogeological data from the coal mine, the average buried depth of the static water level in the mining area was 10 m, and there was a certain hydraulic connection between the two aquifers. Therefore, a water head pressure of 1.8 MPa varying in gradient along the gravity direction was applied at the bottom of the first karst aquifer, and the left and right lower boundaries were set as impervious boundaries. The flow steady was set for seepage calculation. For comparing the influence of water drainage in the karst aquifer on the overlying strata movement and surface subsidence, a reference numerical simulation model of a karst aquifer without water was also established while keeping the other conditions constant, and the monitoring line was set along the strike on the surface with 30 points to analyze the overlying strata movement and surface subsidence. Table 1 illustrates the physical, mechanical and hydraulic parameters of each overlying strata layer. According to the field-monitoring hydrogeological data from the coal mine, the average buried depth of the static water level in the mining area was 10 m, and there was a certain hydraulic connection between the two aquifers. Therefore, a water head pressure of 1.8 MPa varying in gradient along the gravity direction was applied at the bottom of the first karst aquifer, and the left and right lower boundaries were set as impervious boundaries. The flow steady was set for seepage calculation. For comparing the influence of water drainage in the karst aquifer on the overlying strata movement and surface subsidence, a reference numerical simulation model of a karst aquifer without water was also established while keeping the other conditions constant, and the monitoring line was set along the strike on the surface with 30 points to analyze the overlying strata movement and surface subsidence. Table 1 illustrates the physical, mechanical and hydraulic parameters of each overlying strata layer.

Overlying Strata Movement under the Condition of Karst Aquifer without Water
With the advance of the working face, the overlying strata movement and surface subsidence for the karst aquifer without water in the model are illustrated in Figure 3. When the advance of working face was 50 m, the immediate roof was suspended firstly under the action of self-weight stress, and the separation of strata or even collapses could be observed in the middle part of the overlying strata. With the working face advancing to 100 m, the overlying strata began to move, separate and sink due to the varying physical and mechanical properties of each rock layer, and the rock stratum in the middle of the goaf was collapsed and compacted. In addition, there was a large separation in the rock strata above the middle. With the continuous advancing of the working face as shown in Figure 3c-h, the caving strata in the middle of the goaf was continuously compacted, and large separation fractures appeared in the overlying strata near the open cut and the working face. The cantilever structure of the overlying strata was formed at the open cut and working face due to the support of the coal wall. Therefore, a similar triangular area could be observed between the open cut and the working face. Generally, the overlying strata undergo the processes of deformation, separation and collapse with the generation of three zones in the horizontal and vertical directions.

Overlying Strata Movement under the Condition of Karst Aquifer without Water
With the advance of the working face, the overlying strata movement and surface subsidence for the karst aquifer without water in the model are illustrated in Figure 3. When the advance of working face was 50 m, the immediate roof was suspended firstly under the action of self-weight stress, and the separation of strata or even collapses could be observed in the middle part of the overlying strata. With the working face advancing to 100 m, the overlying strata began to move, separate and sink due to the varying physical and mechanical properties of each rock layer, and the rock stratum in the middle of the goaf was collapsed and compacted. In addition, there was a large separation in the rock strata above the middle. With the continuous advancing of the working face as shown in Figure 3c-h, the caving strata in the middle of the goaf was continuously compacted, and large separation fractures appeared in the overlying strata near the open cut and the working face. The cantilever structure of the overlying strata was formed at the open cut and working face due to the support of the coal wall. Therefore, a similar triangular area could be observed between the open cut and the working face. Generally, the overlying strata undergo the processes of deformation, separation and collapse with the generation of three zones in the horizontal and vertical directions.  Table 2 illustrates the front and rear collapse angles of the overlying strata with the different advance distances. It can be observed that the rear collapse angle of the overlying strata stayed constant at 63°, and that the front collapse angle of the overlying strata increased with the increase in advancing distance at first, and then kept stable at 65° after advancing 300 m. In addition, the front collapse angle of the overlying strata was equal to the rear collapse angle when the advance distance of the working face was 200 m.   Table 2 illustrates the front and rear collapse angles of the overlying strata with the different advance distances. It can be observed that the rear collapse angle of the overlying strata stayed constant at 63 • , and that the front collapse angle of the overlying strata increased with the increase in advancing distance at first, and then kept stable at 65 • after advancing 300 m. In addition, the front collapse angle of the overlying strata was equal to the rear collapse angle when the advance distance of the working face was 200 m. Figure 4 illustrates the surface vertical displacement with the increase in advancing distance of the working face for the karst aquifer without water. It can be observed that the surface subsidence caused by coal mining was relatively small (<0.06 m) until the working face advanced to 150 m. When the advancing distance of the working face was 200 m, the surface subsidence in the middle of the goaf was greater than that of both sides, and the growth rate of the maximum surface subsidence increased significantly until the advancing distance of the working face increased to 300 m. Correspondingly, the maximum surface subsidence reached to 0.62 m. And the maximum point of surface subsidence was close to the middle of goaf. In addition, the surface-subsidence curve was symmetrical with the center of the overlying strata in the goaf. Moreover, the surface subsidence above the open cut was always greater than that of above the working face in the whole advancing process of the working face.  Figure 4 illustrates the surface vertical displacement with the increase in advancing distance of the working face for the karst aquifer without water. It can be observed that the surface subsidence caused by coal mining was relatively small (<0.06 m) until the working face advanced to 150 m. When the advancing distance of the working face was 200 m, the surface subsidence in the middle of the goaf was greater than that of both sides, and the growth rate of the maximum surface subsidence increased significantly until the advancing distance of the working face increased to 300 m. Correspondingly, the maximum surface subsidence reached to 0.62 m. And the maximum point of surface subsidence was close to the middle of goaf. In addition, the surface-subsidence curve was symmetrical with the center of the overlying strata in the goaf. Moreover, the surface subsidence above the open cut was always greater than that of above the working face in the whole advancing process of the working face. In terms of the vertical displacement of different overlying strata layers, the curves were generally symmetrically distributed as shown in Figure 5, and the vertical displacement of the overlying strata in the goaf decreased with the increase in the distance away from the coal seam, while the vertical displacement of the overlying strata above the open cut and the coal wall slightly increased with the increase in the distance away from the coal seam. To be specific, the maximum vertical displacement of the overlying strata away from coal seam with 10 m, 34 m and 84 m was 1.79 m, 1.65 m and 0.94 m, respectively.   In terms of the vertical displacement of different overlying strata layers, the curves were generally symmetrically distributed as shown in Figure 5, and the vertical displacement of the overlying strata in the goaf decreased with the increase in the distance away from the coal seam, while the vertical displacement of the overlying strata above the open cut and the coal wall slightly increased with the increase in the distance away from the coal seam. To be specific, the maximum vertical displacement of the overlying strata away from coal seam with 10 m, 34 m and 84 m was 1.79 m, 1.65 m and 0.94 m, respectively.  Figure 4 illustrates the surface vertical displacement with the increase in advancing distance of the working face for the karst aquifer without water. It can be observed that the surface subsidence caused by coal mining was relatively small (<0.06 m) until the working face advanced to 150 m. When the advancing distance of the working face was 200 m, the surface subsidence in the middle of the goaf was greater than that of both sides, and the growth rate of the maximum surface subsidence increased significantly until the advancing distance of the working face increased to 300 m. Correspondingly, the maximum surface subsidence reached to 0.62 m. And the maximum point of surface subsidence was close to the middle of goaf. In addition, the surface-subsidence curve was symmetrical with the center of the overlying strata in the goaf. Moreover, the surface subsidence above the open cut was always greater than that of above the working face in the whole advancing process of the working face. In terms of the vertical displacement of different overlying strata layers, the curves were generally symmetrically distributed as shown in Figure 5, and the vertical displacement of the overlying strata in the goaf decreased with the increase in the distance away from the coal seam, while the vertical displacement of the overlying strata above the open cut and the coal wall slightly increased with the increase in the distance away from the coal seam. To be specific, the maximum vertical displacement of the overlying strata away from coal seam with 10 m, 34 m and 84 m was 1.79 m, 1.65 m and 0.94 m, respectively.   Figure 6 illustrates the overlying strata movement for the karst aquifer with water. The separation and collapse of the overlying strata could be firstly observed in the middle position when the advance distance of the working face was 50 m, and a further large separation between the sandy mudstone and the bottom of the first karst aquifer was revealed due to the influence of mining stress and aquifer seepage when the advance distance of the working face was 100 m. With the advance distance of the working face increasing to 150 m, the separation range between the sandy mudstone and the bottom of the first karst aquifer continued to expand, and then this separation was compacted when the advance distance of working face was 200 m, which could not be observed in the overlying strata movement for the karst aquifer without water. Meanwhile, a large separation between the bottom of the second karst aquifer and the mudstone was generated. With the continuous advancing of the working face, the separated strata were gradually compacted again.

Overlying Strata Movement under the Condition of Karst Aquifer with Water
Mathematics 2022, 10, x FOR PEER REVIEW 8 of 18 separation between the sandy mudstone and the bottom of the first karst aquifer was revealed due to the influence of mining stress and aquifer seepage when the advance distance of the working face was 100 m. With the advance distance of the working face increasing to 150 m, the separation range between the sandy mudstone and the bottom of the first karst aquifer continued to expand, and then this separation was compacted when the advance distance of working face was 200 m, which could not be observed in the overlying strata movement for the karst aquifer without water. Meanwhile, a large separation between the bottom of the second karst aquifer and the mudstone was generated. With the continuous advancing of the working face, the separated strata were gradually compacted again.  Table 3 illustrates the collapse angle of the overlying strata for the karst aquifer with water. Similarly, the rear collapse angle was still 63°, while the front collapse angle increased at first and then kept constant with the increase in the workface advancing, and the front collapse angle was equal to the rear collapse when the advance distance of the working face was 250 m. Generally, the overburden collapse angle in this condition was roughly the same for the karst aquifer without water under each advancing distance of the working face.  Figure 7 illustrates the surface vertical displacement with the increase in advancing distance of the working face for the karst aquifer with water. It can be observed that the surface subsidence caused by coal mining was relatively small (<0.07 m) until the working face advanced to 150 m. When the advancing distance of the working face was 200 m, the surface subsidence in the middle of the goaf was greater than that of both sides, and the  Table 3 illustrates the collapse angle of the overlying strata for the karst aquifer with water. Similarly, the rear collapse angle was still 63 • , while the front collapse angle increased at first and then kept constant with the increase in the workface advancing, and the front collapse angle was equal to the rear collapse when the advance distance of the working face was 250 m. Generally, the overburden collapse angle in this condition was roughly the same for the karst aquifer without water under each advancing distance of the working face.  Figure 7 illustrates the surface vertical displacement with the increase in advancing distance of the working face for the karst aquifer with water. It can be observed that the surface subsidence caused by coal mining was relatively small (<0.07 m) until the working face advanced to 150 m. When the advancing distance of the working face was 200 m, the surface subsidence in the middle of the goaf was greater than that of both sides, and the growth rate of the maximum surface subsidence increased significantly until the advancing Mathematics 2022, 10, 169 9 of 17 distance of the working face increased to 300 m. Correspondingly, the maximum surface subsidence reached 0.72 m, which was greater than that of the karst aquifer without water.
Mathematics 2022, 10, x FOR PEER REVIEW 9 of 18 growth rate of the maximum surface subsidence increased significantly until the advancing distance of the working face increased to 300 m. Correspondingly, the maximum surface subsidence reached 0.72 m, which was greater than that of the karst aquifer without water. In terms of the vertical displacement of the different overlying strata layers for the karst aquifer with water, the curves were generally symmetrically distributed as shown in Figure 8, and the vertical displacement of the overlying strata in the goaf decreased with the increase in the distance away from the coal seam, while the vertical displacement of the overlying strata above the open cut and the coal wall slightly increased with the increase in the distance away from the coal seam. Compared with the condition of the karst aquifer without water, the vertical displacement in the different rock layers for the karst aquifer increased. To be specific, the maximum vertical displacement of the overlying strata away from coal seam with 10 m, 34 m and 84 m was 1.87 m, 1.74 m and 1.2 m, respectively. As shown in Figure 9, the surface subsidence for the karst aquifer with water was greater than the karst aquifer without water. To be specific, the maximum vertical displacement was 0.62 and the subsidence coefficient was 0.305 in the condition of the karst aquifer without water, while the corresponding values were 0.72 m and 0.35, respectively, for the karst aquifer with water. This can be explained by the fact that the water in the karst aquifer flowed to the working face along the fracture with the advance of the working face because the fracture in the rock strata was connected to the karst aquifer. In addition, the drainage of the working face led to a decrease in osmotic pressure of the karst aquifer and the increase in the effective stress of the fractured limestone. In terms of the vertical displacement of the different overlying strata layers for the karst aquifer with water, the curves were generally symmetrically distributed as shown in Figure 8, and the vertical displacement of the overlying strata in the goaf decreased with the increase in the distance away from the coal seam, while the vertical displacement of the overlying strata above the open cut and the coal wall slightly increased with the increase in the distance away from the coal seam. Compared with the condition of the karst aquifer without water, the vertical displacement in the different rock layers for the karst aquifer increased. To be specific, the maximum vertical displacement of the overlying strata away from coal seam with 10 m, 34 m and 84 m was 1.87 m, 1.74 m and 1.2 m, respectively.
Mathematics 2022, 10, x FOR PEER REVIEW 9 of 18 growth rate of the maximum surface subsidence increased significantly until the advancing distance of the working face increased to 300 m. Correspondingly, the maximum surface subsidence reached 0.72 m, which was greater than that of the karst aquifer without water. In terms of the vertical displacement of the different overlying strata layers for the karst aquifer with water, the curves were generally symmetrically distributed as shown in Figure 8, and the vertical displacement of the overlying strata in the goaf decreased with the increase in the distance away from the coal seam, while the vertical displacement of the overlying strata above the open cut and the coal wall slightly increased with the increase in the distance away from the coal seam. Compared with the condition of the karst aquifer without water, the vertical displacement in the different rock layers for the karst aquifer increased. To be specific, the maximum vertical displacement of the overlying strata away from coal seam with 10 m, 34 m and 84 m was 1.87 m, 1.74 m and 1.2 m, respectively. As shown in Figure 9, the surface subsidence for the karst aquifer with water was greater than the karst aquifer without water. To be specific, the maximum vertical displacement was 0.62 and the subsidence coefficient was 0.305 in the condition of the karst aquifer without water, while the corresponding values were 0.72 m and 0.35, respectively, for the karst aquifer with water. This can be explained by the fact that the water in the karst aquifer flowed to the working face along the fracture with the advance of the working face because the fracture in the rock strata was connected to the karst aquifer. In addition, the drainage of the working face led to a decrease in osmotic pressure of the karst aquifer and the increase in the effective stress of the fractured limestone. As shown in Figure 9, the surface subsidence for the karst aquifer with water was greater than the karst aquifer without water. To be specific, the maximum vertical displacement was 0.62 and the subsidence coefficient was 0.305 in the condition of the karst aquifer without water, while the corresponding values were 0.72 m and 0.35, respectively, for the karst aquifer with water. This can be explained by the fact that the water in the karst aquifer flowed to the working face along the fracture with the advance of the working face because the fracture in the rock strata was connected to the karst aquifer. In addition, the drainage of the working face led to a decrease in osmotic pressure of the karst aquifer and the increase in the effective stress of the fractured limestone. Mathematics 2022, 10, x FOR PEER REVIEW 10 of 18 Figure 9. Comparison of surface vertical displacement.

Calculation Model of Surface Subsidence Caused by Aquifer Drainage
When considering the consolidation drainage of a soil aquifer, the compression deformation is usually calculated by formula (1) [28]. The corresponding principle of effective stress can produce an increment of Δp when the pore-water pressure in the soil decreases Δp and the axial deformation of limestone aquifer (εz) is related to various factors such as osmotic pressure drop (Δp), initial seepage pressure (p0), confining pressure (σ3) and other factors in the calculation of drainage deformation of the karst aquifer. However, the osmotic pressure in the fractured limestone can decrease Δp, and the stress increment is different due to the difference in initial osmotic pressure (p0) and confining pressure (σ3), which did not agree with the principle of effective stress. Therefore, the principle effective stress modified by the Biot coefficient αb can be used in the calculation of consolidation-drainage deformation of the karst aquifer, and Equation (1) can be changed to Equation (2).  (2) where αb is the Biot coefficient; εz is the axial displacement deformation; av is the compressibility coefficient of rock and soil mass; e0 is the initial void ratio of rock and soil mass; Δp is the pore pressure.
According to the characteristics of aquifer drainage, the consolidation-drainage-deformation model of the karst aquifer was established as shown in Figure 10.

Calculation Model of Surface Subsidence Caused by Aquifer Drainage
When considering the consolidation drainage of a soil aquifer, the compression deformation is usually calculated by formula (1) [28]. The corresponding principle of effective stress can produce an increment of ∆p when the pore-water pressure in the soil decreases ∆p and the axial deformation of limestone aquifer (ε z ) is related to various factors such as osmotic pressure drop (∆p), initial seepage pressure (p 0 ), confining pressure (σ 3 ) and other factors in the calculation of drainage deformation of the karst aquifer. However, the osmotic pressure in the fractured limestone can decrease ∆p, and the stress increment is different due to the difference in initial osmotic pressure (p 0 ) and confining pressure (σ 3 ), which did not agree with the principle of effective stress. Therefore, the principle effective stress modified by the Biot coefficient α b can be used in the calculation of consolidation-drainage deformation of the karst aquifer, and Equation (1) can be changed to Equation (2).
where α b is the Biot coefficient; ε z is the axial displacement deformation; a v is the compressibility coefficient of rock and soil mass; e 0 is the initial void ratio of rock and soil mass; ∆p is the pore pressure. According to the characteristics of aquifer drainage, the consolidation-drainagedeformation model of the karst aquifer was established as shown in Figure 10.

Calculation Model of Surface Subsidence Caused by Aquifer Drainage
When considering the consolidation drainage of a soil aquifer, the compression deformation is usually calculated by formula (1) [28]. The corresponding principle of effective stress can produce an increment of Δp when the pore-water pressure in the soil decreases Δp and the axial deformation of limestone aquifer (εz) is related to various factors such as osmotic pressure drop (Δp), initial seepage pressure (p0), confining pressure (σ3) and other factors in the calculation of drainage deformation of the karst aquifer. However, the osmotic pressure in the fractured limestone can decrease Δp, and the stress increment is different due to the difference in initial osmotic pressure (p0) and confining pressure (σ3), which did not agree with the principle of effective stress. Therefore, the principle effective stress modified by the Biot coefficient αb can be used in the calculation of consolidation-drainage deformation of the karst aquifer, and Equation (1) can be changed to Equation (2). where αb is the Biot coefficient; εz is the axial displacement deformation; av is the compressibility coefficient of rock and soil mass; e0 is the initial void ratio of rock and soil mass; Δp is the pore pressure.
According to the characteristics of aquifer drainage, the consolidation-drainage-deformation model of the karst aquifer was established as shown in Figure 10. "a" was set as the buried depth of the groundwater level before aquifer drainage. If the rock mass below z = a was saturated, the element stress of plane "dcdh" at depth h can be expressed as follows where γ s is the volume weight of rock mass above the groundwater level; γ sat is the volume weight of saturated rock mass. Pore-water pressure can be expressed as follows where γ w is the volume weight of pore water. The particle stress of rock mass skeleton can be expressed as follows The pore-water pressure caused by the drainage of aquifer will be borne by the skeleton particles of the rock mass. According to the principle of effective stress modified by the Biot coefficient α b , the stress increment of the rock-mass skeleton particles can be expressed as follows According to Equation (2), the small compression deformation of element "dcdh" under the action of stress increment (α b ∆p) can be expressed as follows According to the random medium theory, after the drainage of rock mass, the movement of rock mass above the unit is dW water (x) due to the small subsidence (ds) caused by the drainage of karst aquifer. Meanwhile, the movement of rock mass above the unit is dW(x) + dW water (x) due to the micro subsidence (dh), where dW(x) is the micro unit of surface subsidence caused by the rock mass without considering drainage.
Therefore, the corresponding expression can be expressed as follows According to the random medium theory, the surface forms a small unit subsidence basin. In the drainage water of rock and soil mass, any unit can produce slight compression, and the sinking formula of a water micro unit can be shown as follows where r w is the main influence radius of surface subsidence caused by drainage. According to the principle of superposition, the surface subsidence caused by aquifer drainage can be expressed as follows According to Equations (4)-(8), the formula can be obtained as follows According to the literature, when the drainage of rock mass is not considered and the mining width of the coal seam is "c", the surface subsidence can be expressed as follows [28] where W 0 is the maximum surface subsidence, r is the main influence radius of surface subsidence, m is the mining height, and α is the dip angle of coal seam. Using the probability integral function erf, Equation (12) can be transformed into the following formula.
where erf is the probability integral function, erf The probability integral function can be obtained by looking up the probability integral table.
By substituting Equation (13) into Equation (11), the following formula can be obtained.
It can also be expressed as follows The change of the geotechnical void ratio can be expressed as follows And the surface subsidence caused by underground mining and aquifer drainage can be expressed as follows where W (x) is the surface subsidence caused by underground mining and aquifer drainage, W(x) is the surface subsidence caused by mining obtained by probability integral method.

Theoretical Calculation
According to the above numerical-simulation results, due to the influence of mining, the water flow in the overlying aquifer was recharged downward to the lower aquifer as shown in Figure 11, and there was a hydraulic connection among the three rock strata.
where W0 is the maximum surface subsidence, r is the main influence radius of surface subsidence, m is the mining height, and α is the dip angle of coal seam.
Using the probability integral function erf, Equation (12) can be transformed into the following formula.
where erf is the probability integral function, . The probability integral function can be obtained by looking up the probability integral table.
By substituting Equation (13) into Equation (11), the following formula can be obtained.
[ ] It can also be expressed as follows The change of the geotechnical void ratio can be expressed as follows And the surface subsidence caused by underground mining and aquifer drainage can be expressed as follows where W'(x) is the surface subsidence caused by underground mining and aquifer drainage, W(x) is the surface subsidence caused by mining obtained by probability integral method.

Theoretical Calculation
According to the above numerical-simulation results, due to the influence of mining, the water flow in the overlying aquifer was recharged downward to the lower aquifer as shown in Figure 11, and there was a hydraulic connection among the three rock strata.  However, with the advancement of workface mining, the following processes can change in turn: When the water-conducting fissure zone extended to the first karst aquifer, the first karst aquifer was drained. There was a head difference between the first karst aquifer and the second karst aquifer. When the initial seepage gradient of the rock and soil strata was reached, the overflow from phreatic water to confined water could be observed.
With the expansion of the water flow seepage to the impermeable mudstone, the overall permeability of the impermeable mudstone increased. The second karst aquifer replenished the first karst aquifer through the impermeable mudstone, and the water volume increased. The amount of water supplemented by the second karst aquifer to the first karst aquifer through the impermeable mudstone increased as well.
When the overall water impermeability of the mudstone stratum reduced to a certain value, the first and second karst aquifers could connect with each other and the drawdown of the water level increased with the increase in the thickness of aquifer, which intensifies the surface subsidence value caused by the consolidation-drainage deformation of aquifer. Therefore, the two overlying karst aquifers can be regarded as one aquifer for calculation.
According to field data in the Daojiao coal mine, without considering the drainage of aquifer, the calculation parameters of surface subsidence after coal mining can be expressed as follows.
The coefficient of mining subsidence: q = 0.305. The main influence tangent angle: where ϕ is the internal friction angle (30 • ).
The main influence radius of surface subsidence after C 3 coal-seam mining. × 400) + 1 = 0.13m The maximum value of surface subsidence caused by underground mining and aquifer drainage can be expressed as follows: W max = W water (x) + W 0 = 0.73 m.
According to the calculation results, the surface subsidence caused by aquifer drainage accounted for 17.8% of the total surface subsidence, indicating that the increment of surface subsidence caused by multi-aquifer drainage cannot be ignored.

Observation Scheme
Laser-tracking technology has the advantages of a high measurement accuracy and efficiency, and it is widely used in mine-surface measurement. Therefore, laser-tracking technology was used to observe the surface subsidence with the advance of the 10,301 working faces in the Daojiao coal mine. In order to reduce the impact of temperature fluctuation on the observation results, the observation points arranged on the surface were buried 300 mm below the topsoil. The observation data adopted the national standard of the third-class vertical control point as the reference point in order to ensure the accuracy of the observation results. Figure 12  the third-class vertical control point as the reference point in order to ensure the accuracy of the observation results. Figure 12 illustrates the arrangement of surface-observation stations in the 10,301working faces. An observation line was designed along the trend of the working face with a length of 500 m (baseline A) and an observation point was arranged every 20 m on this line (total 26 observation points). In addition, a total of nine observation lines were arranged along the inclination of the working face with a length of 200 m and observation points arranged every 20 m on these lines. Specifically, observation line B was located above the coal pillar and away from the boundary by 25 m, and observation lines C1 to C8 were located above the goaf and away from open cut by 25

Results and Analysis
The laser-tracking technology was used to observe the surface subsidence during the whole mining period (about one month), and total data records were compiled for each measuring point. Figure 13 illustrates the surface-subsidence curves in the trend and inclination directions of the 10,301 working faces. With the continuous advancement of the working face, the maximum value of the surface subsidence gradually moved to the central position of the goaf, and the surface subsidence increased with the increase in workface advancement until the workface had advanced to 300 m. The maximum surface subsidence was 0.74 m and the subsidence curve was symmetrical with the middle of the goaf as the center. In addition, the surface subsidence near the open cut was always greater than that near the working face. This can be explained by the fact that the affected range of surfaces in the mining process became larger with the advancement of the workface and the cracks caused by mining were gradually compacted due to the collapsed rock filling in the goaf as the roof support. Moreover, the consolidation-drainage of aquifer near the open cut was greater than that near the working face, and the stress originally borne by pore water was transferred to the rock mass to increase the effective stress, which caused the consolidation of rock and the surface subsidence. It can be illustrated that the surface subsidence for the karst aquifer with water was the joint action of coal-seam mining and aquifer drainage.

Results and Analysis
The laser-tracking technology was used to observe the surface subsidence during the whole mining period (about one month), and total data records were compiled for each measuring point. Figure 13 illustrates the surface-subsidence curves in the trend and inclination directions of the 10,301 working faces. With the continuous advancement of the working face, the maximum value of the surface subsidence gradually moved to the central position of the goaf, and the surface subsidence increased with the increase in workface advancement until the workface had advanced to 300 m. The maximum surface subsidence was 0.74 m and the subsidence curve was symmetrical with the middle of the goaf as the center. In addition, the surface subsidence near the open cut was always greater than that near the working face. This can be explained by the fact that the affected range of surfaces in the mining process became larger with the advancement of the workface and the cracks caused by mining were gradually compacted due to the collapsed rock filling in the goaf as the roof support. Moreover, the consolidation-drainage of aquifer near the open cut was greater than that near the working face, and the stress originally borne by pore water was transferred to the rock mass to increase the effective stress, which caused the consolidation of rock and the surface subsidence. It can be illustrated that the surface subsidence for the karst aquifer with water was the joint action of coal-seam mining and aquifer drainage.

Discussion
Based on the numerical simulation and theory analysis, the maximum surface subsidence of the mining area for the karst aquifer without water drainage was 0.62 m and 0.6 m, respectively, and the maximum surface subsidence of the mining area with considering water drainage was 0.72 m and 0.73 m, respectively, which was consistent with the field monitoring result of 0.74 m. It verified the results and the accuracy of the numerical simulation and theoretical calculation. The proposed theory model in this study was

Discussion
Based on the numerical simulation and theory analysis, the maximum surface subsidence of the mining area for the karst aquifer without water drainage was 0.62 m and 0.6 m, respectively, and the maximum surface subsidence of the mining area with considering water drainage was 0.72 m and 0.73 m, respectively, which was consistent with the field monitoring result of 0.74 m. It verified the results and the accuracy of the numerical simulation and theoretical calculation. The proposed theory model in this study was based on the effective stress principle modified by the Biot coefficient α b and can effectively and accurately obtain the surface subsidence while considering water drainage in a karst aquifer and also provide great help in understanding karst topography.
However, there were also some limitations in illustrating the surface subsidence and the overlying strata movement with respect to the water drainage of aquifer. Specifically, the hydrogeological conditions were simplified in the numerical simulation. However, the lithology and hydrogeological conditions of the overlying strata were complex and changeable in the process of real mining environment. In addition, the drainage tests of aquifer rock samples can be further performed to analyze the relationship between axial deformation and the relevant factors (e.g., pore pressure drop ∆p, initial seepage pressure p 0 and confining pressure σ 3 ). Moreover, more than two observation lines can be laid out along the strike of workface to improve the accuracy of the observation results under the condition of economic rationality.

Conclusions
In this study, by taking 10,301 working faces of the Daojiao coal mine in Guizhou Province as the engineering background, comprehensive research methods (e.g., numerical simulation, theory analysis, field monitoring) were adopted in order to illustrate the laws of overlying strata movement and the mechanism of surface subsidence with/without considering water drainage in multi-karst aquifers. The main conclusions are as follows: (1) Compared with the karst aquifers without water, the movement and deformation characteristics of overburden for the karst aquifers with water drainage were quite different, while the collapse angle of overburden was roughly the same. (2) With the advance of the workface, the fracture caused by mining activities can penetrate to a karst aquifer and the water in the aquifer can flow to the workface along the fracture. The water drainage can decrease the osmotic pressure of aquifer, increase the effective stress of fractured limestone and the compression of aquifer, which can intensify the surface subsidence. (3) The maximum surface subsidence of mining area for the karst aquifer without water drainage was about 0.6 m, and the maximum surface subsidence of the mining area when considering water drainage was about 0.73 m. Therefore, the surface subsidence caused by drainage of multi-aquifer accounted for 17.8% of the total surface subsidence. (4) Based on the effective stress principle modified by the Biot coefficient α b , the axial deformation of aquifer with considering water drainage can be obtained, and the field-monitoring results of surface subsidence can also verify the accuracy of the theory-model results.