Numerical Simulation and Application of Transient Electromagnetic Detection Method in Mine Water-Bearing Collapse Column Based on Time-Domain Finite Element Method

: Mine ﬂoods caused by water-bearing collapse columns are one of the most signiﬁcant types of coal mine disasters in advanced detection. The transient electromagnetic method (TEM), one of the most effective means, can detect the collapse column in coal mines. However, the research on TEM used on the ground is mainly limited to half-space that does not follow underground mining conditions. This research used the whole-space TEM to detect the water-bearing collapse columns in front of the roadway. Based on the time-domain ﬁnite element method (TDFEM) and actual coal-measure strata data, this paper built the whole-space geo-electric models, summarized response characteristics, and proposed “geometric gravity center” effect. The results showed that the simulation based on TDFEM could accurately reﬂect the TEM response law, and the advanced detection using TEM in a mine environment was consistent with the actual situation.


Introduction
The development of deep resources has become an unavoidable approach for many nations around the world due to the gradual loss and depletion of superficial mineable resources [1][2][3]. However, as the mining depth of the mine increases, not only is the upper pressurized water level of the stope high, but also the water damage to the mine floor is severe [4]. Different water sources enter the mine, causing water inrush disasters because of the changes in the properties of the rock mass brought on by deep, high in-situ stress, high ground temperature, and high fluid pressure, together with the development of fault fracture zones in the strata themselves and the water-conducting fracture zones generated by mining [5][6][7]. Because of its rapidity, economy and convenience, geophysical methods, widely used to identify hidden water conduction structures such as waterbearing columns, are an effective way to prevent deep well water disasters and solve the problem of pressurized water mining [8][9][10][11]. The transient electromagnetic method (TEM), which has the advantages of high sensitivity to low-resistance bodies and can be easily constructed, is one of the primary methods for detecting water-bearing structures down in the mines [12][13][14][15]. Through the three-dimensional numerical simulation, the study of the transient electromagnetic response of water-bearing columns can not only provide a theoretical foundation for its three-dimensional inversion and data interpretation, but it can also encourage the development and use of this method, which is of great practical significance for the prevention and control of water in coal mines and safe production.
The transient electromagnetic method used in mine hydrogeological exploration is mainly divided into the ground and underground [16][17][18][19]. Affected by the narrow space of the underground tunnel, the detection devices used on the ground are different from those used underground. As shown in Figure 1, the equipment on the floor is a large ring with varying lengths of side from tens of meters to hundreds of meters, while the underground equipment is small with a side length below 3 m. In addition, the geological anomalies to be detected are underground, so the TEM used on the ground has the characteristics of a halfspace response. As shown in Figure 2, the reaction of underground TEM will be affected by the whole space. Unfortunately, most current studies on transient electromagnetic detection are concentrated in the half-space. With the increase of mining depth, the distance between the ground detection device and the detection target will also increase, resulting in a decrease in the resolution of ground detection. Therefore, it is significant to study the whole space as the research object of the underground TEM data processing method.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 2 of 16 widely used to identify hidden water conduction structures such as water-bearing columns, are an effective way to prevent deep well water disasters and solve the problem of pressurized water mining [8][9][10][11]. The transient electromagnetic method (TEM), which has the advantages of high sensitivity to low-resistance bodies and can be easily constructed, is one of the primary methods for detecting water-bearing structures down in the mines [12][13][14][15]. Through the three-dimensional numerical simulation, the study of the transient electromagnetic response of water-bearing columns can not only provide a theoretical foundation for its three-dimensional inversion and data interpretation, but it can also encourage the development and use of this method, which is of great practical significance for the prevention and control of water in coal mines and safe production. The transient electromagnetic method used in mine hydrogeological exploration is mainly divided into the ground and underground [16][17][18][19]. Affected by the narrow space of the underground tunnel, the detection devices used on the ground are different from those used underground. As shown in Figure 1, the equipment on the floor is a large ring with varying lengths of side from tens of meters to hundreds of meters, while the underground equipment is small with a side length below 3 m. In addition, the geological anomalies to be detected are underground, so the TEM used on the ground has the characteristics of a half-space response. As shown in Figure 2, the reaction of underground TEM will be affected by the whole space. Unfortunately, most current studies on transient electromagnetic detection are concentrated in the half-space. With the increase of mining depth, the distance between the ground detection device and the detection target will also increase, resulting in a decrease in the resolution of ground detection. Therefore, it is significant to study the whole space as the research object of the underground TEM data processing method.   widely used to identify hidden water conduction structures such as water-bearing columns, are an effective way to prevent deep well water disasters and solve the problem of pressurized water mining [8][9][10][11]. The transient electromagnetic method (TEM), which has the advantages of high sensitivity to low-resistance bodies and can be easily constructed, is one of the primary methods for detecting water-bearing structures down in the mines [12][13][14][15]. Through the three-dimensional numerical simulation, the study of the transient electromagnetic response of water-bearing columns can not only provide a theoretical foundation for its three-dimensional inversion and data interpretation, but it can also encourage the development and use of this method, which is of great practical significance for the prevention and control of water in coal mines and safe production.
The transient electromagnetic method used in mine hydrogeological exploration is mainly divided into the ground and underground [16][17][18][19]. Affected by the narrow space of the underground tunnel, the detection devices used on the ground are different from those used underground. As shown in Figure 1, the equipment on the floor is a large ring with varying lengths of side from tens of meters to hundreds of meters, while the underground equipment is small with a side length below 3 m. In addition, the geological anomalies to be detected are underground, so the TEM used on the ground has the characteristics of a half-space response. As shown in Figure 2, the reaction of underground TEM will be affected by the whole space. Unfortunately, most current studies on transient electromagnetic detection are concentrated in the half-space. With the increase of mining depth, the distance between the ground detection device and the detection target will also increase, resulting in a decrease in the resolution of ground detection. Therefore, it is significant to study the whole space as the research object of the underground TEM data processing method.   The transient electromagnetic forward method is an effective way to study the transient electromagnetic response, mainly including the time-domain finite-difference (TDFD) method, the finite volume method, the integral equation method, and the finite element method (FEM) [20][21][22][23]. Compared with the above methods, the main advantages of the finite element method are its discrete spatial flexibility and high simulation accuracy for complex geometries [24]. In addition, with the rapid development of computer hardware, the problem of extensive computing time and memory requirements has been effectively alleviated, so the FEM is considered a common method. Kuo and Cho [25] used the FEM to simulate the TEM response of low resistivity ore bodies and discovered that the three-dimensional numerical simulation is essentially comparable with the field measured results, demonstrating the viability of the approach. Since then, Pridmore et al. [26] also carried out exploratory research on the finite element method for solving the forward problem of three-dimensional transient electromagnetic method. Sugeng [27]  on the concept of separating the anomalous field and the background field. Through the TDFEM and tetrahedral mesh, Yin et al. studied the impact of terrain on the transient electromagnetic response [28]. Li et al. used unstructured tetrahedral meshes to study the propagation laws of short electromagnetic fields in 3D forward models of fixed and moving coils [29]. Few studies have been conducted on the numerical simulation of water-bearing collapse columns in three dimensions under full-space conditions, and their response mechanism has not been adequately explained.
Against the above background, the contributions of this study can be summarized as follows: (1) Based on the TDFEM, the whole-space transient electromagnetic responses of different water-bearing collapse columns were simulated by COMSOL. (2) By analyzing multiple sets of numerical simulation results, the "geometric gravity center" effect of the transient electromagnetic response was proposed, which guided the advanced interpretation data of the TEM. (3) The simulation and practical application indicated that the advanced detection could accurately reflect the position of the collapse column by the underground TEM.

Time-Domain Electromagnetic Field Equation
In the application of the TEM to a mine's multi-turn small loops, the electromagnetic field frequency is low, and the displacement current is negligible, in which case the timedomain Maxwell equation of its transient electromagnetic field can be expressed as: where h(r, t) and e(r, t) are the magnetic and electric fields at a certain point r in the whole space. σ(r) and µ(r) are the conductivity and permeability respectively. µ 0 is generally taken as the magnetic conductivity of vacuum. For a multi-turn small loops in mines, By taking the rotation of Equation (2) and substituting it into Equation (1), we can obtain the following expression: Equation (4) is the equation of the magnetic field diffusion.

Vector Finite Cell Method Discretized to the Space Domain
In the calculation of the electromagnetic factors using the three-dimensional finite element method, the calculation area is generally divided into finite tetrahedron or hexahedron elements [30]. The calculation principle of these units can be found in the literature [30]. First, we set up at bedding V e unit in the area, and the four vertices of the tetrahedron are numbered i = 1, 2, 3, and 4, corresponding to the magnetic field strengths on the edge h e i (t). Figure 3 shows the tetrahedral unit structure.
The magnetic field at any point P in the tetrahedral unit E can be expressed as: where N e i (r) is the vector interpolation base function, and h e i (t) represents the magnetic field value on the ith edge in the tetrahedral element E.

Vector Finite Cell Method Discretized to the Space Domain
In the calculation of the electromagnetic factors using the three-dimensional finite element method, the calculation area is generally divided into finite tetrahedron or hexahedron elements [30]. The calculation principle of these units can be found in the literature [30]. First, we set up at bedding unit in the area, and the four vertices of the tetrahedron are numbered i = 1, 2, 3, and 4, corresponding to the magnetic field strengths on the edge ℎ ( ). Figure 3 shows the tetrahedral unit structure. The magnetic field at any point P in the tetrahedral unit E can be expressed as: where ( ) is the vector interpolation base function, and ℎ ( ) represents the magnetic field value on the ith edge in the tetrahedral element E.

Time Discretization
In this paper, the backward Euler method is used to discretize in the time domain using the unconditionally stable backward Euler method, which can be seen as a firstorder backward differential formula. The calculation principle is as follows [31]: In addition, considering the actual situation, the analogue time range is wide, and the current shutdown time is small. In the current shutdown time, the current and magnetic field vary rapidly, and the time step is relatively small. At other times, these parameters are relatively stable, and the time step can be appropriately longer. Therefore, this paper adopts different time steps for different stages, which can help significantly improve the computing efficiency under the premise of ensuring the calculation accuracy.

Boundary Condition Setting
For the calculation of transient electromagnetic fields, the commonly used approximate boundary conditions are the Neumann boundary conditions, Dirichlet boundary conditions, perfect magnetic conductor, absorption boundary conditions, and perfect matching layer (PML). Notably, since Neumann boundary conditions, Dirichlet boundary conditions, and perfect magnetic (electric) conductor all assume that the electromagnetic field value at the boundary nodes is approximately 0, the ideal state of this type of boundary condition is that the calculation region is infinite. For a 3D FEM, this type of boundary

Time Discretization
In this paper, the backward Euler method is used to discretize in the time domain using the unconditionally stable backward Euler method, which can be seen as a first-order backward differential formula. The calculation principle is as follows [31]: In addition, considering the actual situation, the analogue time range is wide, and the current shutdown time is small. In the current shutdown time, the current and magnetic field vary rapidly, and the time step is relatively small. At other times, these parameters are relatively stable, and the time step can be appropriately longer. Therefore, this paper adopts different time steps for different stages, which can help significantly improve the computing efficiency under the premise of ensuring the calculation accuracy.

Boundary Condition Setting
For the calculation of transient electromagnetic fields, the commonly used approximate boundary conditions are the Neumann boundary conditions, Dirichlet boundary conditions, perfect magnetic conductor, absorption boundary conditions, and perfect matching layer (PML). Notably, since Neumann boundary conditions, Dirichlet boundary conditions, and perfect magnetic (electric) conductor all assume that the electromagnetic field value at the boundary nodes is approximately 0, the ideal state of this type of boundary condition is that the calculation region is infinite. For a 3D FEM, this type of boundary condition is often inappropriate because of the constraints of the calculation efficiency and calculation region.
Absorption boundary conditions are the best choice for simulating the propagation of electromagnetic fields in infinite space with limited computational areas and meshes. Because the complete absorption of an electromagnetic field is not realistic, an exact matching layer was proposed by Berenger to absorb electromagnetic waves that propagate outward [32]. Compared to conventional boundary conditions, a PML is equivalent to an area added outside the model that absorbs all the outgoing waves (ideally). This method is perfect for the FDTD truncation of electromagnetic waves; however, because the PML boundary conditions proposed by Berenger do not satisfy Maxwell's equation, it is difficult to discretize finite cells. An infinite cell is a unit that tends to be infinite geometrically and is the perfect complement to the finite cell for solving infinite-domain problems [33]. Therefore, combined with the infinite element method and finite element method, and with the advantages of the PML method, an infinite-element domain is added to the outer boundary of the model, through coordinate mapping, and the magnetic field propagation is simulated to an infinite distance. Because Astley mapping infinite elements can easily handle transient problems, we used the Astley infinite meta-domain method as the boundary condition in this paper [34,35].

Non-Uniform Meshing
In this paper, the computational area is divided using unstructured tetrahedron meshes. The mesh is fine-grained in the emission backline and nearby areas to match the geometry of the area near the small return source and the rapid changes in the electromagnetic field. Due to the difference in the resistivity, the cross-interface of the difference in the semi-spatial resistivity is also refined. In other areas, because of the large geometry and small electrical differences, the tetrahedron mesh gradually expands in the direction of the calculation boundary to ensure the calculation accuracy, which significantly improves the calculation efficiency. Figure 4a,b shows the meshing effect.
the PML boundary conditions proposed by Berenger do not satisfy Maxwell's equation, it is difficult to discretize finite cells. An infinite cell is a unit that tends to be infinite geometrically and is the perfect complement to the finite cell for solving infinite-domain problems [33]. Therefore, combined with the infinite element method and finite element method, and with the advantages of the PML method, an infinite-element domain is added to the outer boundary of the model, through coordinate mapping, and the magnetic field propagation is simulated to an infinite distance. Because Astley mapping infinite elements can easily handle transient problems, we used the Astley infinite meta-domain method as the boundary condition in this paper [34,35].

Non-Uniform Meshing
In this paper, the computational area is divided using unstructured tetrahedron meshes. The mesh is fine-grained in the emission backline and nearby areas to match the geometry of the area near the small return source and the rapid changes in the electromagnetic field. Due to the difference in the resistivity, the cross-interface of the difference in the semi-spatial resistivity is also refined. In other areas, because of the large geometry and small electrical differences, the tetrahedron mesh gradually expands in the direction of the calculation boundary to ensure the calculation accuracy, which significantly improves the calculation efficiency. Figure 4a,b shows the meshing effect.

Verification of Algorithm Accuracy
To verify the accuracy of the finite unit algorithm, the geoelectrical model proposed in this paper is based on a 3D sphere model considering the good symmetry of a sphere.

Verification of Algorithm Accuracy
To verify the accuracy of the finite unit algorithm, the geoelectrical model proposed in this paper is based on a 3D sphere model considering the good symmetry of a sphere. The sphere model has a radius of 180 m, and the 18 m thick spherical shell at its outermost edge is set as the infinite-element domain. The geoelectric model is built inside, as shown in Figure 4a. The procedure is as follows: A circular coil with a radius of 1 m and a crosssectional radius of 25 mm is placed at the center of the sphere, with a coil emission current waveform of 2 A, a shutdown time of 20 µs, a transmit backline of 60, and a receiving line of 40. Using non-uniform meshing, it is divided into 114,834 tetrahedron meshes. The analogue time starts at 0 s, and the current starts to shut down at 0.5 ms, with a shutdown time of 20 µs. The 1.5 ms analogue time ends when the current is switched off. The method of local-refinement time-step sectioning is applied, where a time step of 2 µs is used in the current shutdown process and in the range of 20 s before and after, and a time step of 20 µs is used for the rest of the duration. This can also ensure the calculation accuracy and significantly improve the computational efficiency. Figure 5a,b compares the numerical values obtained from the full-space model and the analytical solution. The compared parameter is the induced electromotive force at the center of the coil. The full-space numerical simulation results agree the analytical solution in each period, and the relative error in the induced electromotive force is below 4.5%.

Whole-Space Transient Electromagnetic Responses of Water-Bearing Collapse Column
First, this paper establishes a sphere space modeled at 180 m and sets the 18 m thick outer boundary layer as an infinite meta-domain. A layered geoelectrical model is established, divided into three layers, representing roof, coal seam, and floor. The specific geometry and geoelectrical parameters are shown in Table 1. For the roadway model of the coal seam, we set the length of the tunnel being dug as 100 m and the roadway section

Whole-Space Transient Electromagnetic Responses of Water-Bearing Collapse Column
First, this paper establishes a sphere space modeled at 180 m and sets the 18 m thick outer boundary layer as an infinite meta-domain. A layered geoelectrical model is established, divided into three layers, representing roof, coal seam, and floor. The specific geometry and geoelectrical parameters are shown in Table 1. For the roadway model of the coal seam, we set the length of the tunnel being dug as 100 m and the roadway section dimensions as 4 m × 4 m. The interior area of the laneway is set as air, the resistivity of which is 1 × 10 5 Ω·m. The peripheral surrounding rock resistivity of the geoelectrical model is set to 100 Ω·m. A wire device serves as the excitation source, and the coil is taken as a circular coil with a radius of 1 m, the transmitting coil is 60 turns, the receiving coil is 40 turns, the coil is placed on the tunnel head, the detecting direction is in front of the tunnel head. The coil and the tunnel head have a 25 mm interval to be consistent with the actual situation. An abnormal geologic body is placed 20 m from the tunnel head, with a resistivity of 1 Ω·m; its parameters are listed in Table 2. To study the transient electromagnetic response law of the collapse column, this paper establishes a collapse column model with an upper surface radius of 10 m, a radius of 30 m on the lower surface, a height of 100 m, a resistivity of 1 Ω·m, and placed 20 m from the excavation head. Figure 6a,b shows the specific model structure. Figure 7a-f shows the inductive electromagnetic field in the collapse column at different times after the current is switched off. When there is no current, the far-left part of the middle of the collapse column is the first to sense the disappearance of the field and produces an inductive electromagnetic field in the area. At this point, the inductive electromagnetic field peak is approximately 4.5 × 10 −9 T, and in 10 µs, the response of the electromagnetic field spreads to the inside of the collapse column and continues to spread to the center area. At 200 µs, the response spreads to the center area of the collapse column, after which it stabilizes, with only the attenuation of the electromagnetic field strength.

Responses of Different Inclination Collapse Column
Previous research has shown that the transient electromagnetic response is first generated by the point sensing of the abnormal body closest to the coil after being shut down and gradually migrates internally over time and eventually to the central part. Therefore,

Responses of Different Inclination Collapse Column
Previous research has shown that the transient electromagnetic response is first generated by the point sensing of the abnormal body closest to the coil after being shut down and gradually migrates internally over time and eventually to the central part. Therefore, this paper proposes a "geometric gravity center" effect of the transient electromagnetic response, i.e., regardless of the morphological distribution of the low-impedance bodies, the place of the first inductive electromagnetic field is the part closest to the coil, and the final area is always concentrated near its geometric gravity center.
To verify the "geometric gravity center" effect, this section established the collapse column models with different inclination angles, including −90, −60, −30, 0, 30, 60, and 90 • . In an actual geological environment, the development of the vast majority of collapse columns is vertical or similar to a vertical formation distribution. The upper surface radius of the collapse column model in each figure is 10 m. The lower surface radius is 20 m, the height is 80 m, and the apparent resistivity is 1 Ω·m. Figure 8a-g shows the magnetic induction strength distribution of the collapse column at different angles on the xoz plane 10 µs after shutdown. After the shutdown, the area closest to the coil of the collapse column is the first to perceive the disappearance of a field, and an inductive electromagnetic field is generated in this area. The collapse column inclination is different, and the position and magnetic induction intensity of the first transient electromagnetic response is different. After a shutdown for 500 µs (Figure 9a-g), the transient electromagnetic response law in the collapse column at various angles is stabilized. column models with different inclination angles, including −90, −60, −30, 0, 30, 60, and 90°. In an actual geological environment, the development of the vast majority of collapse columns is vertical or similar to a vertical formation distribution. The upper surface radius of the collapse column model in each figure is 10 m. The lower surface radius is 20 m, the height is 80 m, and the apparent resistivity is 1 Ω·m. Figure 8a-g shows the magnetic induction strength distribution of the collapse column at different angles on the xoz plane 10 μs after shutdown. After the shutdown, the area closest to the coil of the collapse column is the first to perceive the disappearance of a field, and an inductive electromagnetic field is generated in this area. The collapse column inclination is different, and the position and magnetic induction intensity of the first transient electromagnetic response is different. After a shutdown for 500 μs (Figure 9a-g), the transient electromagnetic response law in the collapse column at various angles is stabilized.   Figure 10a shows that the decay rate of the induced electromotive force of each model, from large to small, follows the order: −90° > −60° > −30° > 0° > 30° > 60° > 90°. The early and late induced electro-

Practical Application
This part seeks to use the above-mentioned response law to direct the field detection and data interpretation of 3114 working face in Mataihao coal mine, thereby verifying the law's accuracy. According to the three-dimensional seismic data of the mine, the DF25 fault run through 3112, 3114, and 3116 working faces, the transport roadway of 3114 working face, and the return airway of 3116 working surface crossed the fault. To find out the distribution of water-rich area of DF25 fault near the transport roadway crossheading of 3114 working face and the return airway crossheading of 3116 working face, it is necessary to carry out underground TEM advanced detection of fault water-richness.
The design height of the 3114 working face is 6 m, 3 1 coal is mined, the buried depth of the coal seam is about 387 m, the thickness of the coal seam is 6~7.7 m, and the average thickness is 6.53 m. The thickness of the coal seam roof is 164.07~255.13 m for Jurassic strata, and the average thickness is 187.47 m. The 3114 face coal seam has a crack-to-production ratio of around 27 times, meaning that when the working face is mined to a height of 6 m, the development height of the water-conducting fracture zone is 162 m. After the recovery of the 3114 working face, the water-conducting fracture zone developed in the Jurassic stability formation formation of the coal seam roof, and was not connected with the Cretaceous Zhidan group sandstone aquifer. The direct water filling source during the mining process of the 3114 working face is the sandstone fracture aquifer of the Yan'an Formation and Zhiluo Formation on the roof of the coal seam, the aquifer is weak in water richness, the normal water inrush of the working face is 100 m³/h, and the maximum inrush is 200 m³/h. The working surface is generally filled with low water intensity. Figure 11 is the result of advanced detection before grouting. In the pictures, the bluegreen and yellow-red regions represent the low resistivity anomaly areas (relatively rich

Practical Application
This part seeks to use the above-mentioned response law to direct the field detection and data interpretation of 3114 working face in Mataihao coal mine, thereby verifying the law's accuracy. According to the three-dimensional seismic data of the mine, the DF25 fault run through 3112, 3114, and 3116 working faces, the transport roadway of 3114 working face, and the return airway of 3116 working surface crossed the fault. To find out the distribution of water-rich area of DF25 fault near the transport roadway crossheading of 3114 working face and the return airway crossheading of 3116 working face, it is necessary to carry out underground TEM advanced detection of fault water-richness.
The design height of the 3114 working face is 6 m, 3 1 coal is mined, the buried depth of the coal seam is about 387 m, the thickness of the coal seam is 6~7.7 m, and the average thickness is 6.53 m. The thickness of the coal seam roof is 164.07~255.13 m for Jurassic strata, and the average thickness is 187.47 m. The 3114 face coal seam has a crack-to-production ratio of around 27 times, meaning that when the working face is mined to a height of 6 m, the development height of the water-conducting fracture zone is 162 m. After the recovery of the 3114 working face, the water-conducting fracture zone developed in the Jurassic stability formation formation of the coal seam roof, and was not connected with the Cretaceous Zhidan group sandstone aquifer. The direct water filling source during the mining process of the 3114 working face is the sandstone fracture aquifer of the Yan'an Formation and Zhiluo Formation on the roof of the coal seam, the aquifer is weak in water richness, the normal water inrush of the working face is 100 m 3 /h, and the maximum inrush is 200 m 3 /h. The working surface is generally filled with low water intensity. Figure 11 is the result of advanced detection before grouting. In the pictures, the blue-green and yellow-red regions represent the low resistivity anomaly areas (relatively rich water area) and the relatively high resistivity areas (weak water area or no water area), respectively. In the direction of roof and layered detection, the abnormal area is located on the left side of the roadway 40 m~right side 40 m, detection direction 30~50 m. After the exploration, the mine conducted drilling design and construction combined with geophysical exploration results. During the drilling construction, the water yield of No. Z2-8 borehole in 3116 return airway was the largest, reaching 66 m 3 /h, and the water yield of other boreholes was 1-11 m 3 /h. The drilling water output area is in the low resistivity anomaly area delineated by mine transient electromagnetic.
the exploration, the mine conducted drilling design and construction combined with geophysical exploration results. During the drilling construction, the water yield of No. Z2-8 borehole in 3116 return airway was the largest, reaching 66 m 3 /h, and the water yield of other boreholes was 1-11 m 3 /h. The drilling water output area is in the low resistivity anomaly area delineated by mine transient electromagnetic.
After determining the specific location of the water-bearing collapse column, multiple boreholes are arranged on the ground for grouting treatment to prevent flooding. When the grouting treatment is completed, the mine TEM was used again for detection. The detection results are shown in Figure 11. As shown in Figure 11, when there is a waterconducting structure, the resistivity value at the anomaly is reduced, and the contour distribution is distorted, deformed, or densely striped. Figure 12 showed that the isoline of apparent resistivity changes stably and uniformly without significant decrease. This showed that after grouting treatment, the collapse column and surrounding rock strata are fully integrated into the relevant electrical characteristics.  After determining the specific location of the water-bearing collapse column, multiple boreholes are arranged on the ground for grouting treatment to prevent flooding. When the grouting treatment is completed, the mine TEM was used again for detection. The detection results are shown in Figure 11. As shown in Figure 11, when there is a waterconducting structure, the resistivity value at the anomaly is reduced, and the contour distribution is distorted, deformed, or densely striped. Figure 12 showed that the isoline of apparent resistivity changes stably and uniformly without significant decrease. This showed that after grouting treatment, the collapse column and surrounding rock strata are fully integrated into the relevant electrical characteristics.

Discussion
Figure 7a-f show the induced electromagnetic field in the collapse column at different times after breaking the current, which is consistent with the research law obtained by Chang's [36] research: After the current is turned off, the induced electromagnetic field in the collapse column first appears at the position closest to the transmitting coil, and then gradually diffuses to the central area and stabilizes. By comparing Figure 8a-g in Section 3.2, each model of the magnetic induction intensity peak from large to small is: −90° = 90° > 60° > −60° > 30° > −30° > 0. When an electromagnetic induction field is generated in the collapse column after the shutdown, the minimum distance between the collapse column and the coil is different due to the difference in the collapse column inclination. The smaller the minimum distance, the greater the magnetic induction strength. In addition, comparing Figure 9a-g, it can be found that the transient electromagnetic responses of all models are finally concentrated near the geometric gravity center of the collapse column. Finally, the law of numerical simulation may lead the data processing and interpretation of the coal mine transient electromagnetic method, as shown by the real field application. The aforementioned response laws can be confirmed by other engineering practice applications in the future, and new transient electromagnetic response laws under complex settings can be researched based on TDFEM.

Discussion
Figure 7a-f show the induced electromagnetic field in the collapse column at different times after breaking the current, which is consistent with the research law obtained by Chang's [36] research: After the current is turned off, the induced electromagnetic field in the collapse column first appears at the position closest to the transmitting coil, and then gradually diffuses to the central area and stabilizes. By comparing Figure 8a-g in Section 3.2, each model of the magnetic induction intensity peak from large to small is: −90 • = 90 • > 60 • > −60 • > 30 • > −30 • > 0. When an electromagnetic induction field is generated in the collapse column after the shutdown, the minimum distance between the collapse column and the coil is different due to the difference in the collapse column inclination. The smaller the minimum distance, the greater the magnetic induction strength. In addition, comparing Figure 9a-g, it can be found that the transient electromagnetic responses of all models are finally concentrated near the geometric gravity center of the collapse column. Finally, the law of numerical simulation may lead the data processing and interpretation of the coal mine transient electromagnetic method, as shown by the real field application. The aforementioned response laws can be confirmed by other engineering practice applications in the future, and new transient electromagnetic response laws under complex settings can be researched based on TDFEM.

Conclusions
(1) Based on the TDFEM, it is discovered that the response law of the water-bearing collapse column under full-space conditions is as follows: After the current is switched off, the induced electromagnetic field in the collapse columns first always appear at the point nearest to the transmitting coil, and then progressively diffuses to the middle area and stabilizes; The smaller the minimum distance, the greater the magnetic induction strength.